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abstract 

A hierarchical flux-based finite element method is developed for both 
one- and two-dimensional thermal-structural analyses. Derivation of the finite 
element equations is presented. The resulting finite element matrices 
associated with the flux-based formulation are evaluated in closed-form. The 
hierarchical finite elements include additional degrees of freedom in the 
approximation of the element variable distributions by the use of nodeless 
variables. The nodeless variables offer increased solution accuracy without 
the need for defining actual nodes and rediscretizing the finite element model. 
Thermal and structural responses obtained using the hierarchical flux-based 
method are compared with results obtained from a conventional linear finite 
element method and exact solutions. Results show that the hierarchical flux- 
based method can provide improved thermal and structural solution accuracy 
with fewer elements when compared to results for the conventional linear 

element method. 
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Chapter 1 
INTRODUCTION 

1.1 Overview 

The finite element method provides a valuable technique for structural 
analysis and design. The method is well suited for the analysis of structures 
with various geometries, loadings, and boundary conditions. Although other 
computational techniques are available, the finite element method is usually 
best suited for problems having complex geometries. A thorough evaluation of 
the structural response induced by aerodynamic loading is an important factor 
in the design of aircraft structures. Thermal and structural finite element 
analyses are often required in the design of high-speed aerospace vehicles to 
prevent structural failure and enhance structural performance. 

For high-speed aircraft, severe aerodynamic heating may occur in local 
areas on the body of the vehicle. Nonuniform heating may produce intense 
local thermal gradients. Since thermal stresses are sensitive to thermal 
gradients, a detailed thermal analysis is required to predict accurate 
temperature distributions needed for evaluation of the thermal stresses. The 
finite element model generally needs to be discretized several times to assure 
convergence and accuracy of the thermal and structural solutions. The process 
of discretizing the finite element model can be time consuming for complicated 


structures and can result in an increasing number of degrees of freedom which 
increases the computational expense. An additional time consuming process 
can be incurred in the transfer of data from the thermal analysis to a form 
suitable for input into the structural analysis. Since the finite element method is 
a widely accepted analysis technique, difficulties and inadequacies in applying 
the method have inspired research for improving the accuracy and efficiency of 
the method. 


1 .2 Literature Review 

The finite element method was first introduced in 1956 as a means for 
analyzing complex aircraft structures [1J. Since its inception, the finite element 
method has become one of the most prominent numerical methods for structural 
analysis. More recently, the finite element method has gained wider 
acceptance for the analysis of thermal and fluid problems. The conventional 
formulation of the finite element equations in all three disciplines and the most 
commonly used element interpolation functions for defining the element 
distribution of the unknown dependant variables can be found in reference 2. 

In general, the accuracy of the finite element solution is improved by 
refining the finite element model using consecutively smaller elements until 
there is convergence of the solution. The method for improving solution 
accuracy by decreasing the element size is known as the h-method. A 
commonly used alternative approach, the p-method, redefines the element 
interpolation functions using more nodes with higher-order interpolation 
functions until the solution converges. An integrated thermal-structural finite 
element approach was introduced by Dechaumphai and Thornton [3-5] which 
improves the solution accuracy and computational efficiency for predicting 
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thermal stresses. The integrated thermal-structural finite element method uses 
a nodeless variable formulation, where additional unknown variables are 
included in the assumed element distribution. These nodeless variables are 
associated with quadratic interpolation functions which produce more accurate 
solutions than the conventional linear element formulation. The nodeiess 
variable formulation provides more accurate transient temperature distributions 
by increasing the degrees of freedom of the element without defining additional 
element nodes, and consequently may yield more accurate thermal stress 
predictions without the need for rediscretizing the finite element model. The use 
of nodeless variables can also be referred to as a hierarchical methodology, 
since the formulation reduces to the conventional linear element formulation 
when the nodeless variables are constrained to zero or eliminated. 

Other approaches for improving the finite element method include the 
development of efficient algorithms for generating the finite element equations. 
A Taylor-Galerkin algorithm, first developed by Donea [6-7] for convective 
transport problems, was applied for the analysis of high-speed flows [8-11]. The 
desire for a single methodology to analyze combined, fluid, thermal, and 
structural interactions led to the extension of the Taylor-Galerkin algorithm for 
the thermal and structural finite element formulations [12-13]. An integrated 
fluid-thermal-structural analysis method [14] was developed for the two- 
dimensional analysis of high-speed flow over leading edges of aerospace 
vehicles. A key feature of the Taylor-Galerkin algorithm is the use of the flux- 
based formulation, where the distribution of the flux of the dependent variable is 
assumed in the same form as the distribution of the dependant variable. The 
flux-based formulation leads to finite element matrices that can be evaluated in 
closed form, whereas the conventional finite element formulation requires 
numerical integration. Another benefit of the algorithm is that nonlinear material 
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properties can be included directly and do not require regeneration of the finite 
element matrices. Also, nonlinear boundary conditions can be incorporated 
easily into the analysis algorithm. These benefits of the flux-based 
methodology led to the further extension of the algorithm for the three- 
dimensional thermal-structural analysis of high-speed wing leading edge 
designs [15-16]. Additionally, a standard two-dimensional eight-node higher- 
order element was incorporated with the flux-based finite element method for 
transient thermal analyses [17]. The use of such a higher-order element 
requires defining additional nodes within the element and consequently 
redefining the finite element model. 

1 .3 Objective 

The objective of this thesis is to develop an improved finite element 
method for predicting accurate thermal and structural responses of structures. 
As an alternative to using higher-order elements, which requires redefining the 
finite element model with additional nodes, this thesis develops and 
investigates the use of nodeless variable finite elements with the flux-based 
finite element formulation. The hierarchical flux-based finite elements have the 
potential of offering a more efficient means of obtaining an accurate thermal- 
structural solution. Both the one- and two-dimensional hierarchical flux-based 
elements are developed for thermal and structural analyses. Transient thermal 
and quasi-static structural analysis capabilities have been developed and are 
contained in a common computer program. The thermal finite element model 
and its temperature solution are completely compatible with the structural finite 
element model. No transfer or manipulation of data is required to obtain the 
thermal loading used in the structural analysis. The finite element results are 
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compared with results obtained using EAL (Engineering Analysis Language 
[18]), a general purpose finite element code frequently used for the thermal and 

structural analyses of aircraft structures. 

Details of the flux-based finite element method for thermal and structural 
analyses are presented in Chapter 2. The basic concepts are introduced along 
with the benefits of the algorithm as compared to the conventional finite element 
formulation. The concept of nodeless variable finite elements, which were 
developed in reference 3 using the conventional finite element formulation, is 
introduced in Chapter 3. Later in this chapter, a one-dimensional hierarchical 
thermal-structural flux-based finite element analysis method is developed and 
results of the method are presented. The methodology is extended to two- 
dimensional elements in Chapter 4 with the development of a membrane 
thermal-structural analysis capability. A summary of the results and concluding 
remarks concerning the effect on accuracy in using the hierarchical flux-based 
finite element method for thermal and structural analyses is presented in 
Chapter 5. 
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Chapter 2 

FLUX-BASED FINITE ELEMENT FORMULATION 
2.1 Basic Concepts 

A common approach to the formulation of many thermal-structural 
problems is to assume that the thermal and structural analyses are uncoupled, 
and that the structural analysis is quasi-static. The thermal and structural 
analyses are uncoupled by neglecting the mechanical deformation rate that 
could alter the temperature in the heat transfer equation. In addition, by 
neglecting the inertia term in the structural equation of motion, the structural 
analysis is deemed quasi-static. These assumptions are credible when the 
temperature change is slow and the coupling effect is negligible. The 
assumptions allow the transient thermal analysis to be performed first, and the 
series of resulting temperatures are then used in performing the structural 
analysis. This approach was applied in the development of the Taylor-Galerkin 
algorithm with a flux-based formulation for thermal-structural analyses [12]. 
The flux-based formulation allows the finite element matrices to be evaluated in 
closed form, and distinguishes the method from the conventional formulation, 
where numerical integration is required. In addition, the flux-based formulation 
has desirable attributes that make it effective for analyzing large transient 
thermal-structural problems, where the thermal analysis can be nonlinear. First, 
the time required to form the finite element matrices using the closed form 
expressions is considerably less than the conventional method, which requires 
numerical integration [15]. Additionally, nonlinear effects, such as temperature 
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dependent material properties, are easily included in the analysis, and do not 
require regeneration of the finite element matrices as in the conventional 
formulation. Because these advantages improve computational efficiency over 
conventional methods, the flux-based algorithm is considered to be more 
suitable for solving complex problems. The development of the finite element 
formulation for transient thermal and quasi-static structural analysis using the 
flux-based formulation is presented in this chapter as a prelude to the extension 
of the flux-based formulation to nodeless variable elements in the following 
chapters. A comparison to the conventional method is made for the two- 
dimensional thermal analysis formulations to present the effect of the flux-based 
formulation on the finite element equations. 

2.2 Governing Equations 


2.2.1 Thermal Analysis 

For an uncoupled thermal-structural analysis, the energy equation 
describing the transient thermal response of the structure in two dimensions can 
be written in the form, 


<HJj dEj DFj 
9t + dx + dy 


= Hj 


( 2 . 1 ) 


where Ut is related to the internal energy, the subscript T denotes the thermal 
analysis, Ej and Fj are the heat transfer components in the x- and y-coordinate 
directions, respectively, and Hj is the internal heat generation. The variable Ut 
and the heat transfer components can be expressed in terms of the temperature 

as, 
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9Ut dT 
dt " pc dt 


E T - Px = -k 0 X (2.2) 

c , dT 

F T = q y - ' k 3y 

where p is the density, c is the specific heat, k is the thermal conductivity, and T 
is the temperature. As shown in equation (2.2), the components of heat flux in 
the two coordinate directions (q x and q y ) are assumed to be related to 
temperature gradients by Fourier's law. Both the specific heat and the thermal 
conductivity may be temperature dependent. 


2.2.2 Structural Analysis 

The quasi-static structural response is governed by the equilibrium 
equations that can be written in the form, 


5{E S } + d{F s } 
ax + dy =0 


(2-3) 


where the subscript S denotes the structural analysis. For two-dimensional 
problems the vectors {Es} and {Fs} contain the stress components given by 

{Es} T = [ o x x X y ] 

{Fs} T = [ x X y Oy ] (2 4) 
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where the stress components o x , o y , and x X y are assumed to be related to 
displacement gradients and the temperature by generalized Hooke's law, 

{o} = [C] {€) + {p} (T-T 0 ) (2-5) 

where {a} is the vector of stress components, {e} is the vector of strain 
components, [C] is the matrix of material elastic constants, {p} is the vector of 
thermal expansion parameters, T is the temperature, and T 0 is the reference 
temperature for a zero stress state. 

2.3 Solution Procedure 


For simplicity and to illustrate the generality of the flux-based algorithm, the 
thermal or structural governing equation is written in the form of a single 
equation as 


au aE aF _ 

ar + S? + a7'° 


( 2 . 6 ) 


where the first term is zero for the quasi-static structural formulation. 

The basic objectives of the Taylor-Galerkin algorithm are to: (1) use a 
Taylor series expansion of the U variable to establish recursion relations, and 
(2) use the Galerkin method of weighted residuals [2] for spatial discretization to 
derive the finite element equations. A more detailed description of the Taylor- 
Galerkin algorithm is presented in the following chapters for the development of 
the nodeless variable flux-based finite element equations. 
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For the thermal analysis formulation, the temperature, T, is the dependent 
variable and is contained in U, where 

au AU AT 

at - At =pc ^T ( 2 - 7 ) 

The term AT is the change in temperature from the previous time step, n, to the 
current time step, n+1, where AT = (T"+i - T n ). The next step in the formulation 
of the finite element equations is to assume a spatial variation of the dependent 
variable throughout the selected element type. A natural coordinate system, 
which simplifies the element geometry, is used to define the spatial variations in 
non-dimensional form. The distribution of temperature, T, for the four-node 
bilinear quadrilateral element, is assumed in the form 

4 

T(x,y,t) = T Nj(C,n) Tj(t) = {N(C,rt)}T {T(t)} (2.8) 

where {N(£,t|)} t is the row matrix of the nodal interpolation functions in natural 
coordinates, £ and q, and {T} is the vector of nodal temperatures. A 
conventional four-node quadrilateral element shape is shown in figure 1(a) and 
the transformation to the natural coordinates is shown in figure 1(b). Since finite 
element matrices are in the form of integrals over element areas, transformation 
to natural coordinates permits the integrals to be evaluated over a square 
region. The Cartesian coordinates are related to the natural coordinates by 


x « [N(Cn>] {x} 

y = [N(C,n>] {y} 


(2.9) 




where {x} and {y} are the vectors of nodal Cartesian coordinates for the element. 
Since the same interpolation functions, Nj, are used to interpolate the 
temperature and the spatial coordinates, the element is called an isoparametric 
element. For the four-node quadrilateral element, the interpolation functions 
are defined by 


Ni=*(K) (1-rj) 
N 2 = |(1+C) (1-ri) 
N 3 = ^(1+C) (1+rt) 
N 4 = |(K)(1+ri) 


The variable, U, which is directly related to T, is also assumed in the same 
form as equation (2.8). As a consequence of the Taylor-Galerkin algorithm, U 
becomes the unknown variable. The thermal finite element equations are 

solved for U at each time step. The temperatures are then determined from 
equation (2.7). 

A feature of the flux-based algorithm, which differs from the conventional 
finite element method, is that the flux-based algorithm expresses the variation of 
the element fluxes E and F in the same form as the element dependent variable 
[6-7], that is 


U(x.y.t) = [N(C,n)]{U(l» 

E(x,y,t) - [N(U)]{E(t» (2.11) 

F(x.y.t) = [N(C,n)]{F(t)) 
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where {U} is the vector of unknown nodal quantities of the variable U. In the 
thermal context, {E} and {F} are the vectors of the element nodal heat fluxes and 
are related to the nodal temperatures by equation (2.2) and equation (2.7). In 
the structural context, {E} and {F} are vectors of the element nodal stresses, 
which are related to the unknown nodal displacements by use of equation (2.5) 

and the strain-displacement relations. 

Application of the Taylor-Galerkin algorithm to the energy equation, 
equation (2.1), and the flux-based assumption [12] results in the transient 
thermal finite element equation 

[M]{AU} n = {R x )i +{Ry)i + {fty? (2-12) 


where [M] denotes the mass matrix 

[M] = 



(2.13) 


which may be diagonalized to produce a lumped mass matrix. The vector {AU} n 

is the change in the nodal values of the variable U between the time steps t n+1 
and t n where t n = nAt. The vectors {Rx}" and {Ry>? are associated with the fluxes 

within the element in the x- and y-coordinate directions and {R}£ is the boundary 

term associated with the flux across the element boundary. These vectors are 
defined by 


{Rx)i = At [D x ] {E} n 


(2.14a) 



where 


( D *I - J^jWdA 

A 

(Ry}?= At(D y ] {F} n 

where 

! D y) - J^ 1 (NpdA 
{Rfa = -At [B] {q} 

where 

[B] = J{N}{N} T dS 


(2.14b) 

(2.15a) 

(2.15b) 

(2.16a) 

(2.16b) 


The vector {q} contains the components of the nodal heat flux normal to the 
element surface boundary. The thermal boundary conditions are applied with 
the vector (q) expressed in equation (2.16a), where the surface nodal heat flux 
q is replaced by the quantities representing any one of several different types of 
boundary conditions 



0 (insulated) 

q s (specified heating) 
h (T-T r ) (surface convection) 

4 4 

eo(T-TJ (surface radiation) 


(2.17) 
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For the structural analysis, the Taylor-Galerkin flux-based formulation 
produces matrices identical to equations (2.14-2.16), where {E} and {F} 
represent the nodal stress components given in equation (2.4). More details of 
the formulation and boundary conditions for the structural analysis are given in 
reference 15. 


2.4 Closed-Form Finite Element Matrices 

All the element integrals obtained from the flux-based formulation can be 
evaluated in closed form. The closed-form matrices apply to the quadrilateral 
as well as the hexahedral element shapes. The availability of closed-form 
expressions eliminates the need for numerical integration in the evaluation of 
element matrices. The closed-form expressions for the flux-based finite element 
matrices, [M], [D x ], [D y ], and [B] defined in equations (2.13-2.16), are merely a 
function of the element geometry. The use of the symbolic manipulation 
program, MACSYMA[19], simplified the evaluation of the closed-form 
expressions. The closed-form expressions for the coefficients in the finite 
element matrix [D x ], for the four-node quadrilateral element are 


Dx (1,1) = * D x (3,3) = - (y 4 - y 2 ) / 6 
D x (2,2) = - D x (4,4) = - ( Y1 - y 3 ) / 6 
D x (1,3) = - D x (3,1) = -(y 4 -y 2 )/12 
D x (2,4) = - D x (4,2) = - (yi - y 3 ) / 12 
D x (1 ,2) = - (y 4 + y 3 - 2y 2 ) / 1 2 
D x (1 ,4) = - (2y 4 + y 3 - y 2 ) / 1 2 

D x (2,1) = (y 4 + y 3 - 2yi) / 12 (2.18) 

D x (2,3) = - (y 4 - 2y 3 + yi) / 12 

1 5 



D x (3,2) = (y 4 - 2y 3 - Yi) / 12 
□x (3,4) = (2y 4 - y 2 - y 1 ) / 12 
Dx (4,1) = -(y 3 + y2-2yi)/12 
D x (4,3) = - (2y 3 - y 2 - yi) / 12 


where Xj and yi are the nodal coordinates. Expressions for [D y ], [M] and [B] are 
similar and are given in reference 15. The flux-based finite element matrices 
are also independent of material properties. This feature removes the necessity 
of reforming the element matrices at every time step when material properties 
are temperature dependent. 

The formulation of the finite element equations using the conventional 
finite element method also begins with the governing equations for heat transfer 
and structural equilibrium expressed in equations (2.1- 2.4). For the thermal 
quadrilateral element formulation, the temperature is also assumed to be in the 
same form as equation (2.8) using the element interpolation functions defined in 
equation (2.10) and the natural coordinate transformation expressed in 
equation (2.9). The thermal finite element equations are in the form 

(Cpl + [K X ](T) + [K y ] (T) = {H} (2.19) 


where {T} is the vector of unknown nodal temperatures, [C p ] is referred to as the 
capacitance matrix, {H} is the internal heat generation vector, and [K x ] and [K y ] 
are the conductance matrices associated with heat conduction in the x- and y- 
coordinate directions, respectively. The matrices are defined by 

(Cp] = f pc {N} {N} T dA (2.20) 

16 


k 


( 2 . 21 ) 


( 2 22 > 

where p is the density, c is the specific heat and k is the thermal conductivity. It 
can be observed from equations (2.20-2.22) that the conventional finite element 
matrices are dependent on the material properties which, in general, are a 
function of the temperature. In addition, there is no closed-form expression for 
[K x ] and [K y ] for arbitrary quadrilateral element shapes. Numerical integration is 
thus required to compute these element matrices. The Gauss integration 
method is used where the integral is expressed as a sum of the weighted terms 
evaluated at Gauss points. Gauss weighting factors and integration points can 
be found in reference 2. In performing numerical integration, the accuracy of 
the matrices depends on the number of Gauss points used. Two Gauss points 
in each coordinate direction are normally used for the bilinear quadrilateral 
element. Due to their dependency on material properties and the need for 
numerical integration, the process of generating the conventional finite element 
matrices can be computationally expensive. 

To predict temperatures and temperature gradients accurately in a 
structure subjected to aerodynamic heating, refined finite element mesh sizes 
may be needed at some locations. Finer meshes, and hence smaller elements, 
require small time steps for analysis solution stability, such as for explicit 
solution algorithms. Hence, the finite element equation needs to be solved 
many times when performing a transient thermal analysis. Temperature 
dependent material properties are also often necessary to model thermal effects 
accurately in a transient thermal analysis. These requirements may make the 
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conventional finite element method computationally expensive for modelling 
large-scale aircraft structures. The flux-based finite element method offers the 
advantages of closed-form finite element matrices and ease in representing 
temperature dependent material properties while providing equivalent solution 
accuracy. Hence, the flux-based finite element method may be more suitable 
for analyzing large nonlinear, transient thermal-structural problems. The 
benefits of the flux-based finite element method inspired the extension to 
nodeless variable elements developed herein. 



Chapter 3 

ONE-DIMENSIONAL NODELESS VARIABLE FINITE ELEMENTS 
USING FLUX-BASED FORMULATION 

3.1 Element Interpolation Functions 

The fundamental basis of the finite element method is that a continuum of 
arbitrary shape can be modeled by an assemblage of simple shapes. For one- 
dimensional problems, the elements are line segments. The line segments are 
assembled to model a one-dimensional domain, which may be either a one- 
dimensional slab or rod continuum. The variation of a dependent variable 
within the element is then approximated as a function of the nodal variables and 
interpolation functions. The conventional finite element method normally 
utilizes linear elements where a linear distribution of a dependent variable 
within an element is assumed using linear interpolation functions. 

For hierarchical finite elements, additional unknown variables, 
sometimes called nodeless variables, are introduced in the assumed 
distribution of a dependent variable for an element to provide a nonlinear 
distribution of the dependent variable. For thermal problems, temperature is the 
dependent variable. A two-node one-dimensional element and typical element 
temperature distributions for the conventional and hierarchical finite elements 
are shown in figure 2. The hierarchical finite element with one nodeless 
variable assumes the element temperature distribution in the form 


2 

» 


1 

Cl 

I— ►x 

K L H 

(a) Two-node one-dimensional element 



(b) Assumed temperature variation in a 
conventional linear element 



(c) Assumed temperature variation in a 
hierarchical nodeless variable element 


Figure 2. One-dimensional thermal finite element and typical 
element temperature distributions. 
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(3.1) 


3 

T(x) = X N i(*) T ' = ( N ( X )) T ro 
i =1 

where {T} is the vector of unknown variables and {N(x)} T is the row vector of 
element interpolation functions. The nodal temperatures are Ti and T 2 , and T 3 
is the nodeless variable. Similarly, in the structural analysis, the element 
displacement, u, is the unknown variable and is expressed in the form 

3 

u(x) = X Ni( x ) u ' = (N(x)} T {u} (3.2) 

i =1 

where {u} is the vector of unknown variables and {N(x)} T is the same row vector 
used to approximate the element temperature distribution. Once again, ui and 
U 2 are nodal displacements, and U 3 is the nodeless variable. The element 
interpolation functions for an element of length L are defined by 

Ni = 1 -[ 


N 2 = l 0.3) 

N 3 = t( 1 -t) 


where Ni and N 2 are the nodal interpolation functions and N 3 is the nodeless 
variable interpolation function. Note that the nodeless variable does not 
represent the actual nodal temperature or displacement, but rather is directly 
related to the magnitude of the nonlinear variation of the element temperature 
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and displacement distributions. When the nodeless variable is constrained to 
be zero, the approximations of the element temperature and displacement 
distributions reduce to the conventional linear element approximation. For the 
hierarchical finite element formulation with one nodeless variable, a quadratic 
distribution of the unknown variables is assumed which is capable of 
representing the general solution more realistically. 

3.2 Derivation of Flux-based Finite Element Equations 

3.2.1 Thermal Analyst 

The transient thermal response of a structure is governed by the energy 
equation. For a one-dimensional, uncoupled, thermal-structural analysis with 
no internal heat generation, the energy equation can be written in the form 

dU dE 

at + dx = 0 (3.4) 

The variable U, and the heat flux, E, are defined by 

ay aT 
at _pc at 

p_ ^ 

fc “ q x = - k ax (3.5) 


where T is the temperature, p is the material density, c is the specific heat, and k 

is the thermal conductivity. Both specific heat and thermal conductivity may be 
temperature dependent. 
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The Taylor-Galerkin algorithm is applied to the governing equation, 
equation (3.4). A Taylor series expansion of the variable U(x,t) in time is 
needed to establish recursion relations. The Taylor series expansion of U(x,t) in 
time is in the form of an infinite series 

u(x,t n+1 ) = u(x,t"> ♦ ™ (t n+1 -t") + ^rfr (t" +1 - 1") 2 

+ 1.|U (t n, , . tn) 3 t ... (3.6) 

The change in the variable U with respect to time is defined as 

AU = U n+1 - U n (3.7) 


For the first order accurate approximation in time, the change in the variable U 
is approximated as 


dU n 

All = At 


(3.8) 


Substituting for the first derivative of the variable U from equation (3.4), equation 
(3.8) becomes 


AU + Atf r = 0 (3.9) 

The Galerkin method of weighted residuals is now applied to minimize the error 
of the approximation of the dependent variable over the element length, 
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(3.10) 


| Nj R dx = 0 


where Nj, the interpolation functions, are used as the weighting functions and R 
is the residual. For the nodeless variable formulation, i = 1 to 3, establishing 
three equations for minimizing the error of the three unknown variables. By 
substituting equation (3.9) for the residual, where R equals the left hand side of 
equation (3.9), equation (3.10) becomes 

J N, AU dx + At /N|^dx=0 ( 311) 

Integration by parts on the second term in equation (3.1 1) yields 

Nj AU dx = At J ^ E" dx - At ( Nj(0) E"(0) - Nj(L) E n (L) ) (3. 1 2) 


The finite element approximations are now needed for discretization in space. 

Since the variable U is directly related to temperature, it follows from equation 
(3.1) that 


where 


3 


AU(x) = £ Nj(x) AUj = {N(x)}T {AU} 
i=1 


AUj = pcAT i = pc(T j n+1 -ff) 


(3.13) 


(3.14) 


The flux-based assumption discretizes the heat flux in the same form as the 
variable U, where heat flux E n replaces AU in equation (3.13). The nodal 
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values of heat flux, Ei and E2, are related to the temperature gradients through 
Fourier's law. Since heat flux is related to the gradient of the temperature, and 
the temperature distribution is assumed to be quadratic, the assumed heat flux 
must be a linear distribution, requiring that the nodeless variable heat flux, E3, 
must be zero. Thus the heat flux distribution reduces to a linear distribution in 
the form 


E"(x) = X Ni(x) E" = {N(x)}T{E}" 
i=1 


(3.15) 


where {E} n is the nodal heat flux vector and { N(x) } T is the row matrix of the 
linear interpolation functions Ni and N2. The vector {AU} and heat flux nodal 
vectors are independent of the integration over the element length. The finite 
element approximations, equations (3.13 and 3.15), are substituted into 
equation (3.12) to yield the finite element equation in the form 

(M]{AU) = At [D] {E} n - At {By} (3. 1 6) 

where [M] denotes the matrix 

[M] = J{N}{N} T dx (3.17) 

The first term of the right hand side of equation (3.16) is associated with the heat 
flux within the element, and {By} is the boundary term associated with the heat 
flux across the element boundaries. The matrix [D] and vector {By} are defined 
by 
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(3.18) 


I 

I 


Pi' jlr 1 < n<*) ) T d* 

L 


f-q(x = 0)' 
{Bt} = | q(x = L) - 

I 0 , 


(3.19) 


where q in equation (3.19) represents the thermal boundary conditions and is 
replaced by any one of the several different types of boundary conditions given 
in equation (2.17). The closed-form expressions for the terms in matrix [MJ and 
[D] are given in Appendix A. 

The nodal heat flux vector in equation (3.16) is related to the nodal 
temperatures through Fourier's law and the finite element approximation for 
temperature. The terms in the vector are evaluated at the corresponding nodes 
where x = 0 at node 1 and x = L at node 2. The vector is defined by 




d{N£ 

dx 


(T} n )node 1 

(T)")„ 0 <le 2 J 


(3.20) 


where | n ode i symbolizes the evaluation of the quantity at node i. Since {E} n is 
dependent on nodal temperatures and the nodal values for thermal 
conductivity, k, it needs to be updated at every time step for transient thermal 
analyses. The matrix [D] is independent of thermal properties and needs only to 
be evaluated once, prior to the transient analysis. The finite element equation, 
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equation (3.16), is solved for AU at every time step. The nodal temperatures are 
determined from equation (3.14). 

3.2.2 Structural Analysis 

The one-dimensional quasi-static structural response is governed by the 
equilibrium equation written in the form 


where E represents the element stress in the x-direction. The element stress is 
defined by 


E = ci - 02 (3.22) 

where oi and 02 are the one-dimensional components of the stress vector 
associated with the displacement gradients and the temperature, respectively. 
The element stress, E, is related to displacement gradients and temperatures by 
Hooke’s law for thermal stress problems. For one-dimensional problems, the 
stress components reduce to 

°<= e ! < 3 - 23 > 

c 2 = E a ( T(x) - T 0 ) (3.24) 

where E is the modulus of elasticity, u is the displacement in the x-direction, a is 
the coefficient of thermal expansion, T(x) is the temperature distribution, and T 0 
is the reference temperature for a zero stress state. 
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As with the thermal formulation, the Galerkin method of weighted 
residuals, equation (3.10) is applied to minimize the error of the approximation 
of the dependent variable over the element length. Here, the governing 
equation, equation (3.21), represents the residual, where R equals the left hand 
side of equation (3.21). Integration by parts is performed to produce an element 
integral term and a boundary integral term for application of applied stress 
boundary conditions. The resulting equation is in the form 

J to Ol d * - f 57 °2 dx - ( Ni(0) 0(0) - Nj(L) o(L) ) = 0 (3.25) 

where i = 1 to 3 yields three equations for minimizing the error of the 
hierarchical finite element approximations. The finite element approximations 
are needed to discretize equation (3.25) in space. The finite element 
approximation of the displacement dependent variable is given in equation 
(3.2). For the structural formulation, the flux-based assumptions are applied to 
the stress components which are analogous to heat flux in the thermal 
formulation. The flux-based assumptions discretize the stress components in 
the same form as the dependent variable. The oi component is a function of the 
displacement gradient, where the displacement distribution is assumed to be 
quadratic. The a i stress component reduces to a linear approximation in the 
same manner that the flux approximation, which is a fuction of the temperature 
gradient, reduces to a linear approximation. Since the 02 stress component is 
directly dependent on temperatures and the temperature distribution 
approximation is quadratic, the approximation for 02 is also quadratic. The 
resulting stress component approximations are thus given by 
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(3.26) 


2 

oi(x) = X Nj(x)jOij = { N(x) } T {ai} 
i=1 

3 

a 2 (x) = X N i( x )°2i = (N(x)} T { 02 } (3.27) 

i=1 

where { N(x) } T is the row matrix of the linear interpolation function, {oi} is the 

vector of nodal values for the stress component associated with the 
displacement gradient, {N(x)} T is the row matrix of the interpolation functions 
given in equation (3.3), and {o 2 } is a vector of known values associated with 

temperature and will be defined subsequently in the thesis. The flux-based 
assumptions, equations (3.26 and 3.27), are substituted into equation (3.25) to 
yield the finite element equation. 

For a one-dimensional, quasi-static, structural analysis, the finite element 
equation is in the form 

[D]{o 1 } = [D 2 ]{o 2 }+{B} (3.28) 

where the matrix [D] is identical to the matrix [D] produced in the thermal 
formulation, and is defined by equation (3.18). The matrix [D 2 ] associated with 
the o 2 stress component and the vector {B} associated with the boundary 

conditions are defined by 

[D 2 ] = jfWcix (3.29) 

L 


29 



-o(x = 0) 
{B} = * c(x = L) 
0 


(3.30) 


where a in equation (3.30) represents an applied stress boundary condition. 

The expressions for the terms in the finite element matrix [D 2 ] are given in 
Appendix A. 

The terms in the vector (o,} are related to nodal displacements through 

equation (3.23) and the finite element approximation tor displacements, 
equation (3.2). The vector (ai) is defined by 



where the gradients of the interpolation functions given in equation (3.3), 
evaluated at the nodes, are substituted into equation (3.31). The displacements 
in equation (3.31) are extracted to yield the {ai} stress component in the form 

(°i) = f p ]{u} (3.32) 

The {o!} nodal stress component vector is given here in terms of the unknown 
displacement vector, {u}, and a matrix, [P], which is a function of the element 
length and the modulus of elasticity. The coefficients in the matrix [P] are given 
in Appendix A. The expression for the stress vector {a,}, is substituted into the 

finite element equation, equation (3.28), which can now be written in terms of 
the unknown displacement vector. 
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The vector {o 2 } is derived by equating the flux-based assumption for o 2 , 
given in equation (3.27), with the definition of o 2 given in equation (3.24), where 
T(x) is replaced by the finite element approximation for T(x), equation (3.1). 
The expression is written in the form 

3 3 

Y Ni(x)o 2i = E a ( £ Ni(x)T, - T 0 ) (3.33) 

i=1 i=1 


Evaluating equation (3.33) at the element nodes where i = 1 and i = 2 , and x 
0 and x = L, respectively, gives the values of the vector {o 2 } as 


(o 2 } = 


Ea (T 1 - To) 1 
< Ea (T 2 - T 0 ) > 
k Ea T3 . 


(3.34) 


The stress component vector, |o 2 ), is a vector of known quantities because the 
element nodal temperatures T, and T 2 . and the element thermal nodeless 
variable T3 are obtained from the thermal analysis. 


3.3 Thermal-Structural Analysis Algorithm 

The hierarchical flux-based finite element analysis method is 
implemented in a computer program for performing transient thermal and quasi- 
static structural analyses. A flow chart for the one-dimensional thermal- 
structural finite element analysis method is shown in figure 3. Once the finite 
element model is input, the finite element matrices [D] and [D 2 ] are assembled 
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Input Finite Element Model I 

r 

Assemble Matrices 

- element [D] and [D 2 I 

- system [M] 

— i 

Thermal Analysis 

1 . compute element vectors {E 0 } and {By} 

2. compute element {RHS} = At([D]{E} n - {By}) 

3. assemble system {RHS} 

4. solve for AU ( [M] {AU} = {RHS} ) 

5. calculate {T} n+1 = {AU/pc 0 } + {T} n 

6. update T n =T n+1 , t n =t n + At 


Structural Analysis 

1 . compute element [LHS] = [D][P] 

2. compute element vectors { 02 } and {B} 

3. compute element {RHS} = [D 2 K 02 } + {B} 

4. assemble system [LHS] and {RHS} 

5. solve for {u} ( [LHS]{u} = {RHS}) 



Figure 3. Flow chart for one-dimensional hierarhical flux-based 
thermal-structural analysis. 








along with the system [M] matrix. The thermal analysis is performed first and 
consists of: (1) computing the element vectors {E n } and {By}; (2) computing the 
element right hand side vectors, where {RHS} is the right hand side of equation 
(3.16); (3) assembling the system {RHS}; (4) solving the system equations for 
{AU}; (5) computing the {T} n+1 vector from equation (3.14); and (6) updating the 
temperatures and time step for proceeding in a transient thermal analysis. If the 
time t n = t s , where t s in the time set for a quasi-static structural analysis, the 
structural analysis begins and the vector {T} n is used for thermal loading. The 
structural analysis consists of: (1 ) computing the element [LHS] matrix, where 
[LHS] represents the left hand side matrix of the finite element equation ([LHS] = 
[D][P]); (2) computing the element vectors { 02 } and {B>; (3) computing the 
element {RHS}, where {RHS} represents the right hand side of equation (3.28); 
(4) assembling the system [LHS] and {RHS}; and (5) solving for the vector {u}. If 
t s = tend, where t en d is the time set to end the thermal-structural analysis, the 
analysis is complete. 


3.4 Applications of One-Dimensional Methodologies 

To evaluate the hierarchical flux-based finite element method, four one- 
dimensional thermal and structural problems are presented. The first two 
example problems are for the transient thermal analysis of a copper slab with 
constant material properties and with temperature dependent material 
properties. The following two example problems are for the thermal and 
structural analysis of a copper rod with one end constrained and with both ends 
constrained. Results obtained by the hierarchical flux-based method are 
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compared with solutions obtained using the conventional method and available 
exact solutions. 


3 - 4 - 1 Iransient Thermal An alysis of a Conner Slah 

A copper slab 0.02 in. in length subjected to an applied heating rate of 
200 Btu/in. 2 -sec on the face at x = 0 and a prescribed temperature equivalent to 
the initial temperature at x = L is analyzed for two cases: (1) constant material 
properties and an initial temperature of Tj = 0 °F; and (2) temperature 
dependent thermal conductivity and an initial temperature of Tj = -410 °F. A 

schematic diagram of the finite element model and the material properties for 
copper are shown in figure 4. 

For case one, with constant material properties, the exact solution to the 
governing equation for one-dimensional transient heat transfer can be derived 
from the method of separation of variables and is in the form of an infinite series 
[20], The exact solution for the transient temperature distribution is given by 


oo 

T(x.»= 9(k*l . M Vi-1)" (gn-fl)it(L-x) ,-k (2n+1) * „ . 

' ' k k „2 Z^(2n+1)2 sm 2L 8 *P (- [ — 21 " 1 0 

n=0 F 1 * 


(3.35) 


where the boundary conditions and initial condition are 


T(L,t) = 0 ( 3 .36) 

T(x,0) = 0 
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^ — L = 0.02 in. — H 


(a) Schematic diagram of thermal 
finite element model 



Temperature, °F 


(b) Material properties for copper 

Figure 4. Schematic diagram of thermal finite element model 
of a copper slab and material properties for copper. 
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The temperature distributions within the slab obtained from the exact solution 
are shown in figure 5 at two times during the transient thermal response. Also 
shown in figure 5 are the finite element solutions obtained from the hierarchical 
flux-based and conventional methods. The flux-based method yields an 
accurate prediction of the exact solution using two nodeless variable elements, 
whereas the conventional method required the use of ten elements to obtain a 
temperature distribution prediction within one percent of the exact solution, not 
shown. The overall average error between the finite element and exact 
temperature distributions is defined by 

Hell D/ 

error = TPflT x 100 % ( 3 -37) 


where 


and 


II e H 2 = l j (Texact *Tfe ) 2 dx 
0 


II Til - ^ J^(Texact) 2 dx 


where T ex act is the exact temperature distribution and Tfe is the finite element 
temperature distribution over the length, L, of the finite element model. Using 
two elements with the conventional method results in a 29% overall average 
error at t = 0.0001 sec. and a 3% overall average error at t=0.001 sec. As 
shown in figure 5, the hierarchical finite element method provides a more 
accurate prediction of the temperature distribution than the conventional 
method for the same number of elements, especially at early times in the 
transient solution when there is a large temperature gradient within the slab. 
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Exact 

2 nodeless variable 



Figures. Case one: Temperature distributions in a copper slab 
with constant material properties. 
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For case two with temperature-dependent thermal conductivity, an initial 
temperature of Tj = -410 <>F is used. At this low initial temperature, the thermal 
conductivity of copper is highly nonlinear as shown in figure 4. The transient 
thermal response of the copper slab, using the temperature dependent 
conductivity shown in figure 4, is predicted using both the hierarchical flux- 
based and conventional methods. Two temperature distributions are shown in 
figure 6 for t = 0.0001 sec. during the transient thermal analysis. Ten 
conventional elements are necessary for the temperature distribution to 
converge to within one percent of the temperature distribution obtained using 
eight conventional elements. The ten-conventional-element temperature 
distribution was then used as a reference solution. The hierarchical flux-based 
method predicted an accurate temperature distribution, within one percent of 
the reference solution, using only four nodeless variable elements. From the 
two examples presented above, the hierarchical flux-based method provides 
accurate temperature distributions using fewer elements. 

3 - 4 * 2 Thermal-Structural Ana lysis of a Popper Rr>ri 

A copper rod one inch in length is analyzed for the transient thermal and 
quasi-static structural response at t = 0.5 sec. A schematic diagram of the 
thermal and structural finite element models for the two cases analyzed are 
shown in figure 7. In both cases, shown in figure 7(a), the thermal analysis 
consists of an applied aerodynamic heating rate of q = 10 Btu/in. 2 -sec at x = 0 
and a prescribed temperature equal to the initial temperature of 70 °F at x = L. 
The temperature distributions obtained using the exact solution, the hierarchical 
flux base method, and the conventional method are shown in figure 8 for t = 0.5 
sec. As can be observed from figure 8, two hierarchical flux-based elements 
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reference solution 

(10 conventional elements) 

4 nodeless variable 



Figure 6. Case two: Temperature distributions in a copper slab 

with temperature dependent conductivity at t = 0.0001 sec. 
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(b) Case one: structural boundary conditions, free end 



(c) Case two: structural boundary conditions, fixed ends 


Figure 7. Schematic diagrams of thermal and structural finite element 
models of a copper rod. 
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Exact 

2 nodeless variable 



Figure 8. Temperature distributions in a copper rod at t = 0.5 sec. 



are sufficient to predict the temperature distribution accurately, whereas four 
conventional elements are required to closely represent the exact solution. 

Two quasi-static structural analysis cases are investigated using the 
temperature distributions shown in figure 8 at t = 0.5 sec. The structural 
boundary conditions for case one, shown in figure 7(b), constrains the end at x 
= 0 and allows for free expansion at x = L. The initial temperature of Tj = 70 °F 
from the thermal analysis was used as the reference temperature for a zero 
stress state, T 0 = 70 °F. The displacement distributions obtained for case one 
using the exact solution, the hierarchical flux-based method, and the 
conventional method are shown in figure 9. Once again, four conventional 
elements are needed to closely approximate the exact solution whereas only 
two hierarchical flux-based elements are sufficient to provide accurate results. 
The next case investigated, case two, assumes both ends are constrained as 

shown in figure 7(c) and a reference temperature of T 0 = 70 °F for a zero stress 
state. 

The displacement distributions and the average element stresses 
obtained for case two are shown in figure 10(a) and figure 10(b), respectively. 
As shown in figure 10(a), two hierarchical flux-based elements were needed to 
closely represent the exact solution, whereas two conventional elements were 
insufficient to predict the displacement distribution accurately. The use of two 
conventional elements underpredicts the maximum displacement by 13 % 
whereas two nodeless variable elements overpredicts the maximum 
displacement by only 3.4 %. Using two conventional elements resulted in a 3% 
error in predicting the stress distribution in the copper rod whereas two 
hierarchical elements predicted the stress distribution within a 1 % error margin. 
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Exact 

2 nodeless variable 

flux-based elements 
— O — 2 conventional elements 
□ 4 conventional elements 



Figure 9. Case one: Displacement distributions in a copper rod 

subjected to thermal loading and constrained at x = 0 (t = 0.5 sec.). 
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Exact 

2 nodeless variable 
flux-based elements 
- - O - — 2 conventional elements 


0 4 conventional elements 




Figure 10. Case two: Structural response for a copper rod with 
both ends constrained. 
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The previous examples demonstrate the ability of the hierarchical method to 
predict accurate thermal and structural responses using fewer elements than 
required by the conventional method to obtain the same accuracy. The 
following chapter includes the extension of the hierarchical flux-based method 
for two-dimensional thermal-structural analyses. 


45 



Chapter 4 

TWO-DIMENSIONAL NODELESS VARIABLE FINITE ELEMENTS 
USING FLUX-BASED FORMULATION 

4.1 Element Interpolation Functions 

In this chapter, the flux-based method is extended to develop a two- 
dimensional finite element formulation using nodeless variables. As with the 
conventional two-dimensional bilinear element formulation, a general 
quadrilateral element shape is employed for formulating the nodeless variable 
element interpolation functions. To simplify the element integrations arising in 
the finite element matrices, the transformation from Cartesian to natural 
coordinates shown in figure 1 is utilized. The relation between the two 
coordinate systems given in equation (2.9) is applied in the development of the 
two-dimensional finite element equations with nodeless variables. 

For the hierarchical method, a nonlinear variation of the dependent 
variable is assumed over the element surface. Preserving the four-noded 
quadrilateral element, the nonlinear variation is established by introducing 
additional degrees of freedom in the approximation of the dependent variable. 
For the thermal analysis formulation, the distribution of temperature over the 
element surface is assumed in the form 
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(4.1) 


8 

T = X Nj(C,Tl) Ti(t) = {N(C,n)) T (T(t)} 

i=1 


where {N(£,ti)} t is the row matrix of element interpolation functions and (T(t)} is 
the vector of unknown variables. The nodal temperatures are Ti through T 4 , 
and the nodeless variables are T 5 through T 8 . For the structural analysis 
formulation, the displacement distributions are assumed in the same form 


8 

u = £ Ni(M Uj = {N(C,ti)}T {u} (4.2) 

i=1 
8 

v = X N i(C.il) v i = {N(C,T1)) T {v} (4.3) 

i=1 


where u and v are the displacements in the x- and y-coordinate directions, 
respectively. The vectors of unknown variables, {u} and {v}, contain the four 
nodal displacements and four nodeless variables. The displacement 
distributions can be expressed in a combined form as 



(4.4) 


where [N] is the combined matrix of interpolation functions for the structural 
formulation and {8} is the vector of nodal displacements and nodeless variables. 
These matrices, [N] and {8}, are defined by 


[N] = 


Ni 0 N 2 0 . . . N 8 0 
0 Ni 0 N 2 . . . 0 N 8 _ 


(4.5) 
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U2 


V 2 


(5} = 


< . 




U8 J 


(4.6) 


The element interpolation functions, N it i = 1 to 4 are identical to those used for 
the conventional bilinear four-node element given in equation (2.10), and Nj, i = 
5 to 8 are the nodeless variable interpolation functions given by 

N s=f (i-C 2 ) (i-n) 

N 6 = |(i+0(i-n 2 ) 

N 7=|(i-C 2 ) (i+n) (47) 

N 8 = |(i-0(i-n 2 ) 


where each interpolation function varies quadratically along one edge and 
vanishes along the other edges of the element. 

Utilizing the nodeless variable interpolation functions in equations (4.1 - 
4.3) provides a quadratic variation of the temperature or displacement 
distribution over the element with only four element nodes. A schematic of 
typical element temperature distributions for the nodeless variable element and 
conventional element are shown in figure 11. The magnitude of the nonlinear 
variation on an element edge depends on the magnitude of the nodeless 
variables. If the nodeless variables are constrained to zero in the 
approximation of the temperature or displacement distribution given in 
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equations (4.1 - 4.3), the distribution reduces to the conventional bilinear 
element approximation. 



T(x, y, t) 



Figure 1 1 . Typical two-dimensional finite element temperature distributions. 
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4.2 Derivation of Flux-based Finite Element Equation 


4.2.1 Thermal Analysis 

The transient thermal response for a two-dimensional uncoupled 
thermal-structural analysis is governed by the energy equation given in 
equation (2.1), where the terms in the equation are defined in equation (2.2). 
The flux-based Taylor-Galerkin algorithm is applied to the governing equation 
to yield the finite element equations. As with the one-dimensional formulation, 
the Taylor series approximation for the variable AU, equations (3.6 - 3.8), is 
utilized to establish recursion relations. The approximation, equation (3.8) is 
substituted in the two-dimensional governing equation, equation (2.1), to yield 

3E n 3F n 

AU+At aT +At = ° (4.8) 


The Galerkin method of weighted residuals is applied to minimize the error of 
the finite element approximations over the element area, 

f Nj R dA = 0 (4.9) 

where Nj, the nodeless variable element interpolation functions, are used as the 
weighting functions and R is the residual. For the quadrilateral nodeless 
variable element, i = 1 to 8, establishing eight equations for minimizing the error 
of the eight unknown variables. The left hand side of equation (4.8) is 
substituted for the residual to yield 

J Nj AU dA + At J Nj ^-dA + At f Nj ^-dA = 0 (410) 
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Integration by parts is performed on the second and third terms in equation 
(4.10) to yield 

J Nj AU dA = At 


^ ^ ^ F 0 dA - At jNj qn dS (4.11) 


where the first two terms on the right hand side represent the heat conduction 
within the element and the last term on the right hand side represents the heat 
flux across the element boundaries. The quantity qn is the heat flux normal to 
the element boundary. 

The next step in the formulation of the thermal finite element equation is 
to discretize equation (4.11) in space by implementing the finite element 
approximations. The variable, AU, being directly related to temperature, is 
discretized in the same form as equation (4.1) and is given by 

8 

AU = X N|(C,H) AUi = {N(C,ti)}T{AU} (4.12) 

i=1 

where 

AUj = P C ATj = p c (Tj n+1 - Tf) (4.13) 


The flux-based assumptions discretize the heat flux in the x- and y-coordinate 
directions in the same form as the dependent variable. As with the one- 
dimensional formulation, the variations of heat flux over the element surface 
reduce to bilinear approximations in the form 

4 

E" = £ NiK,n) E, n = { N(C,n) ) T (E") (4.14) 

i=1 
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(4.15) 


4 


F" - I N|(U) F" ■= 

i=1 


{ N(C,n) } T {F n j 


where { N(C/n) } T is the row matrix of linear interpolation functions given by 

equation (2.10), and {E n } and {F n } are the nodal heat flux vectors in the x- and y- 
coordinate directions, respectively. 

The two-dimensional transient thermal finite element equation can be 
written in matrix form as 

[M]{AU} = At [D x ] {E n } + At [D y ] {F”> - At [B] {q} (4. 1 6) 

where [M] denotes the matrix associated with {AU}, [D x ] and (D y J are associated 

with the heat transfer within the element, and [B] is the boundary matrix. These 
matrices are given by 


(MJ = ^f{N}{N}TdA 

(4.17) 

(D *1 - fd? ( N ) T dA 

A 

(4.18) 


(4.19) 

[B] = J{N}{¥}TdS 

(4.20) 
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The transformation from Cartesian to natural coordinates used to evaluate these 
matrices requires relating the gradients of the interpolation functions in both 
coordinate systems. The chain rule is applied to obtain the relationship 


(du; 

ac 

dN 

lanJ 


> J 



= [j] 


dN' 

dx 

dN 

IdyJ 


(4.21) 


where [J] is the Jacobian matrix. Relating x and y to the natural coordinates 
using equation (2.9), [J] is given by 



Ja_Nj 
i=i ac 

J 3Nj 
i=1 


Xi 


Id Nj 

I x yi 

i=i ac 


« a_Ni 
i=i an 


yi 


’ Jll Jl2 

. J21 ^22 


(4.22) 


From equations (4.21 - 4.22), it follows that 


dx 


f 3N ] 

ac 

► . 1 

" Ml 

J22 

- J 12 1 

P N 1 

ac 

dN 

= [J]’ 1 < 

dN 

. -J21 

Jll J 

> 

dN 

^dyJ 








(4.23) 


where [J]' 1 is the inverse and |J| is the determinate of the Jacobian matrix [J]. 
The derivatives of the interpolation function gradients occurring in equations 
(4.18 - 4.19) with respect to the Cartesian coordinates are replaced by the 
corresponding gradients in natural coordinates as 
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Using the relationship dA = |J| d£ dn. the integration of the matrices given in 
equations (4.17 - 4.19) can now be evaluated with respect to natural 
coordinates over the square element area shown in figure 1. All matrices 
arising in the transient thermal finite element equation for the quadrilateral 
element, equation (4.16), can be evaluating in closed form. The evaluation of 
the matrices, equations (4.17 - 4.19), was greatly simplified by using 
Mathematica [21], a general purpose computer software system for performing 
algebraic manipulation. The closed form expressions for the coefficients in the 
two-dimensional transient thermal finite element matrices and the evaluation of 
|J| are given in Appendix B. 

The nodal heat flux vectors, {En} and {F"}, contained in the transient 
thermal finite element equation, equation (4.16), are related to temperature 
gradients through Fourier's law and are given by 



Lode 1 


'node 2 

' (-^on) 

► 

Lode 3 


node 4 

J 


(4.26) 



{F n } = 




(Jfm) 

(■‘ i S C l T "l) 

(Jfm) 


node 1 


node 2 


node 3 


node 4 


(4.27) 


where the gradients of the interpolation functions with respect to the Cartesian 
coordinates are replaced by the expressions in natural coordinates given in 
equations (4.24 - 4.25). These vectors are a function of nodal temperatures and 
hence need to be updated at every time step for the transient thermal analysis. 
For problems involving temperature dependent conductivity, the thermal 
conductivity is also easily updated. The boundary surface nodal heat flux vector 
{q} in equation (4.16) represents the heat flux normal to the element boundary 
and can be replaced by any of the several different types of boundary heating 
conditions given in equation (2.17). The transient thermal finite element 
equation, equation (4.16), is solved for the nodal change in the variable, {AU}, at 
every time step. The temperatures at the new time step, n+1, are then 
determined using equation (4.13). 

4 . 2.2 Structural Analysis 

The two-dimensional quasi-static structural response is governed by the 
equilibrium equations. By neglecting the body forces, these equations are, 

£*3r-o < 428 > 
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where o x and ay represent normal stresses in the x- and y-coordinate directions 
respectively, and t xy represents shear stress. As with the thermal formulation, 

the Galerkin method of weighted residuals, equation (4.9), is applied to each 
equilibrium equation to yield 

/ | Ni <^ + ^> dA ' 0 (4.30) 

+ 1fc ! > dA = 0 (4.31) 

where Nj, i = 1 to 8 are the nodeless variable interpolation functions given in 
equations (2.10) and (4.7). Integration by parts is performed on each term in 
equations (4.30 and 4.31) to generate the element boundary integrals yielding 
the equations 



jNj T x dS 
> 

(4.32) 

fNjTy dS 

(4.33) 


where T x and T y are surface tractions on the element boundaries in the x- and 
y-coordinate directions, respectively. Equations (4.32 - 4.33) are combined to 
yield sixteen equations for evaluation of the sixteen degrees of freedom of the 



displacement components for the two-dimensional nodeless variable structural 
element. The combined representation is written in the form 


j [B S ] T {a} dA = J [Nf {T s } dS 


(4.34) 


where the subscript s denotes the structural form of the matrices. The matrix 
[N] t is the transpose of the combined interpolation function matrix defined in 
equation (4.5). The matrix, [B S ] T , is given by 



0 


0 

8Ni 

dy 



[BslH * 



0 



dNa 

ay 

9Nb 
9x -I 


(4.35) 


and the vectors {o} and {T s } contain the stress components and surface 
tractions, respectively. These matrices are defined by 


foxl 


{G} = 


1 °y 


r = {^i} - M 


LtxyJ 


(4.36) 


(Ts) = 



(4.37) 
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where {0^ and {a 2 } are the two-dimensional components of the stress vector 
associated with the displacement gradients and the temperature, respectively. 
The first component vector is related to the displacement gradient through the 
generalized Hooke’s law in the form 


M = [CJ{£} 


where 


r 


9 u 

dx 






9v 

dy 




9 u 9 v 

vdy + dx) 


( 4 . 38 ) 


[C] is the matrix of material elastic constants and {e} is the vector of strain 

components. The second stress component vector is related to the temperature 
given by 


M = [C] {a} (T(x,y) - T 0 ) (4.39) 

where {a} is the vector of thermal expansion parameters, T(x,y) is the element 
temperature distribution, and T 0 is the reference temperature for a zero stress 
state. The matrix [C] and vector {a} are dependent on whether the problem 

being analyzed assumes a state of plane stress or plane strain. For the plane 
stress problem and an isotropic material, [C] and {a} are defined by 
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(4.40) 


[C] = 


1-v2 


" 1 v 0 


f \ 



a 

v 1 0 

and {a} s' 1 

a * 

1 

X \°* 

O 

o 


.0, 


and for the plane strain problem [C] and {a} are defined by 


[C] = 


(1+v)(1-2v) 


1-v V 

0 


'a(1+v)' 

v 1-v 

0 

and {a} = < 

a(1+v) 


1-2v 


. 0 , 

° 0 

2 -1 



(4.41) 


where E is the modulus of elasticity, v is Poisson’s ratio, and a is the coefficient 
of thermal expansion. Substituting equation (4.36) into equation (4.34) yields 

f [B s ] t {cti} dA + ^ [Bs] T {°2} dA * J [N S ] T {T s } dS (4.42) 

The finite element approximations are needed to discretize equation 
(4.42) in space. For the structural formulation, the flux-based assumptions 
discretize the element stresses in the same form as the displacement 
discretization given in equations (4.2 - 4.3). It follows from the one-dimensional 
formulation that the first stress component, {oi}, reduces to the bilinear 
approximation. The second stress component, {a 2 }. is directly related to 
temperature which is quadratic. Hence, the flux-based assumptions are defined 
by 


{ai} = [Ni] { ai } 


(4.43) 
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where 


Tn>t {oT {"o"}T ” 

[Nil- {0'}T{¥}T{o‘}T 
_{"0}T {'0}T {¥}T_ 
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and 


({ <*1x }i* 1 ,4) 


{ CTl } = 




({ °1y }i=1 ,4) 


.({ <*1xy 


{02} = [N2] { 02 } 


m T {0} T -1 


[N2I = 


{0) T {N}T 

L {0}T {0}T J 


{ <*2 } = 


j ({ G2x }i*1,8) 
] ({ <*2y }i=1,8) 


1 


^ { 0 } J 


(4.44) 


The vectors { 0 } and {0} are null vectors. The transpose of the vectors is given 
by 
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{ 0 } T = [ 0 0 0 0 ] 

{0} T = [0000000 0] 


(4.45) 


For { oi }, i = 1 to 4 yields the values of the first stress components evaluated at 
node i. For {^2 }, i = 1 to 8 yields the values of the second stress components 

with respect to nodes i = 1 to 4 and the nodeless variables, i = 5 to 8. The terms 
in {&\ } are related to displacement gradients using equation (4.38) and the 

finite element approximation of the displacements, equations (4.2 - 4.3). The 
terms in {”02 } are related to the temperatures using equation (4.39) and the 

finite element approximation of the temperature distribution, equation (4.1). 

The two-dimensional hierarchical flux-based finite element equation is 
written in the form 


[Di] { 01 } = [D 2 ] { 02 } - [B] {T s } (4.46) 


where the matrices [Di], [D 2 ], and [B] are defined by 


Pi] * 

f [B S ]T [N,l dA 

(4.47) 

[D 2 ] = 

f [B S ] T [N Z ] dA 

(4.48) 

[B] = 

f [N] t dS 

(4.49) 


The surface traction vector, {T s }, is assumed here to be constant over the 
element surface. Transformation to natural coordinates is required for 


61 



©valuation of these matrices. As with the thermal matrices, terms in the matrix 
[B S ] T defined in equation (4.35) are replaced by the corresponding gradients in 
natural coordinates given in equations (4.24 - 4.25). Also, using the 
relationship dA = Mld^dq allows for evaluation of the finite element matrices for 
the two-dimensional hierarchical flux-based finite element equation. To 

evaluate the finite element equation in terms of the unknown displacements, the 
vector { ai j is expressed in terms of the displacement vector which can be 

defined as 


{ °i ) = [ p ] {5} (4.50) 

The closed-form expressions for the terms in the two-dimensional structural 
finite element matrices, [Di], [D 2 ], [B], and [P] and the vector } are given in 

Appendix B. Replacing { oi } in equation (4.46) by the relationship given in 
equation (4.50), the finite element equation is evaluated in terms of {5}, the 
vector of unknown nodal displacements and nodeless variables. The values of 
{5} are then used to evaluate the element stresses using equation (4.36) and 
equations (4.38 and 4.39). 

4.3 Applications of Two-Dimensional Methodologies 

To evaluate the two-dimensional nodeless variable flux-based method, 
two example problems are presented. Both example problems are for a plane 
stress analysis of a copper plate. The first problem is a plate with a linear 
temperature distribution. An exact solution for the displacement distributions is 
available for this problem, providing a means for verifying the accuracy of the 
nodeless variable flux-based solution. The second problem analyzed is for a 
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plate with a linear distribution of an applied heat flux over one of the plate 
boundaries. Solutions for both of the example problems are compared with 
solutions obtained using the conventional finite element method. 

A schematic of the structural finite element model for the first problem 
analyzed is shown in figure 12(a). The copper plate is 10 in. long in the x- 
direction and 5 in. wide in the y-direction. A reference temperature for zero 
thermal stress, T 0 = 0, was assumed. The plate is free to expand and is 
constrained to prevent rigid body motion. The temperature distribution within 
the plate was assumed to be a linear function of x, where T(x) = lOx, as shown 
in figure 12(b). Because the plate is subjected to a linear temperature 
distribution and is free to expand, there is no thermal stress. The exact 
displacement distributions are given by 

u(x,y) = 5 a (x 2 - y 2 ) (4-51 ) 

v(x,y) = 10axy (4-52) 

where a is the coefficient of thermal expansion for copper. The exact u- 
displacement distributions at y = 0 and at y = 5 in. are plotted in figures 13(a) 
and 13(b), respectively. Also shown in figure 13(a) and 13(b) are the 
displacement distributions obtained from the nodeless variable flux-based and 
conventional methods. As can be seen from figure 13(a) and 13(b), one 
nodeless variable flux-based element yields the exact displacement 
distributions, whereas one conventional element was insufficient to accurately 
describe the u-displacement distributions. As shown in figure 13(a) and 13(b), 
four conventional elements were required to closely represent the exact 

solution. 
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10 in. 

(b) Assumed temperature distribution in structural 
finite element model 


Figure 12. Schematic diagram of structural finite element model and 
assumed temperature distribution. 
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Exact 

1 nodeless variable 
flux-based element 




Figure 1 3. Displacement distributions in a copper plate 
subjected to a linear temperature distribution. 
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The second problem analyzed consisted of thermal and structural analyses of a 
4 in. square copper plate. The schematic diagram of the thermal and structural 
finite element models are shown in figures 14(a) and 14(b), respectively. For 
the thermal analysis, an initial temperature of Tj = 70 <>F is assumed, where the 
boundary conditions consist of: (1) an applied linearly distributed heating rate, 
q, at x = 4 in., as shown in figure 14(a); (2) a prescribed temperature of 70<> F at 
y = 4 in.; (3) insulated at y = 0; and (4) insulated at x = 0. The temperature 
distributions along the boundaries y = 0 and x = 4 in., obtained using the 
nodeless variable flux-based and conventional methods are shown in figure 
15(a) and 15(b) for the time t = 5 sec. The solution for the 16x16 mesh of 
conventional elements is considered to be the reference solution. A 4 x 4 mesh 
of nodeless variable flux-based elements is required to obtain an accurate 
representation of the reference solution. The temperature distributions at t = 5 
sec., obtained using the 16 x 16 mesh of conventional elements and 4x4 mesh 
of nodeless variable elements, are used as the thermal loading in the structural 
analysis. The thermal finite element model discretizations used to obtain the 
temperature distributions were also used as the structural finite element model 
discretizations. A reference temperature for zero thermal stress, T 0 = 70 °F is 
assumed. As shown in figure 14(b), the structural model is constrained from 
displacement in the y-coordinate direction along the boundary y = o, and 
constrained from displacement in the x-coordinate direction at the corner x = y = 

0. The displacement distributions were obtained using the nodeless variable 
flux-based method, with the 4 x 4 finite element model discretization required to 
accurately represent the reference temperature distribution. The displacement 
distributions were also obtained by using the conventional method, with the 16 
x 16 finite element model discretization. The maximum displacements and 
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(b) Schematic diagram of structural model 

Figure 14 . Schematic diagram of thermal and structural finite 
element model of a copper plate. 
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4x4 nodeless variable 
flux-based elements 




Figure 15. Temperature distributions in a copper plate at t = 5 sec. 
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displacement gradients were observed along the boundary at x = 4 in. for the v 
displacement. The v displacements along the boundary at x = 4 in. obtained 
using the nodeless variable flux-based and conventional method are shown in 
figure 16(a). The maximum stress occurred along the boundary at y = 0 and 
normal to the y-coordinate. The average nodal stress distributions along the 
boundary where y = 0 obtained using the flux-based and conventional method 
are show in figure 16(b). Although only 4x4 nodeless variable elements were 
required to represent accurately the reference temperature and displacement 
distributions, a mesh of 8 x 8 nodeless variable elements is required to 
represent accurately the stress distribution. 
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4x4 nodeless variable 

flux-based elements 
—o— 1 6x1 6 conventional elements 



(a) v displacement distributions at x = 4 in. 


8x8 nodeless variable 
flux-based elements 



(b) Normal stress distribution in y-coordinate 
direction at y = 0 


Figure 1 6. Displacement and stress distributions in a copper plate 
at t = 5 sec. 
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Chapter 5 

CONCLUDING REMARKS 


A hierarchical finite element method using a flux-based formulation 
technique is developed for both one-dimensional and two-dimensional thermal- 
structural analyses. The derivation of the finite element equations is presented 
along with the resulting finite element matrices. The finite element matrices 
associated with the flux-based method can be evaluated in closed-form which 
distinguishes the flux-based method from the conventional finite element 
formulation which requires numerical integration for evaluation of the finite 
element matrices. The hierarchical element is established by introducing 
additional degrees of freedom into the assumed distribution of the unknown 
variables by the use of nodeless variables. Employing hierarchical finite 
elements provides improved solution accuracy without reconstructing new finite 
element models. The technique also allows a single finite element model to be 
used for both the thermal and structural analyses, thus eliminating the difficulty 
in transferring data between the analyses. 

Several thermal and structural example problems are analyzed to 
investigate the ability of the hierarchical flux-based method for predicting 
accurate thermal and structural responses. The resulting solutions obtained by 
using the hierarchical flux-based method are compared with solutions obtained 
using the conventional finite element method and the exact solution when 
available. From the resulting transient thermal solutions, the hierarchical flux- 
based method demonstrates an ability to provide more accurate temperature 
distribution results using fewer elements than the conventional finite element 
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method, especially when large temperature gradients are present in the 
structure. With the quasi-static assumption, the temperature distribution results 
at a specified time can be used as thermal loading in a structural analysis. To 
implement a temperature solution, the corresponding structural model needs to 
have the same discretization as the thermal model. Hence, the conventional 
method usually requires a structural model with more elements then the 
hierarchical flux-based method requires to provide an accurate temperature 
distribution. The hierarchical flux-based finite element method also produced 
accurate structural displacements and stresses using fewer elements than 
required by the conventional method. In general, the hierarchical flux-based 
elements show some improvements in accuracy for predicting thermal and 
structural responses as compared to using an equivalent number of 
conventional elements. 

The hierarchical flux-based method is presently devoloped for the 
analysis of one- and two-dimensional membrane structures. The method is 
general and could be extended to three-dimensional thermal and structural 
analysis capabilities and developed for plate bending problems. 
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Appendix A 

Closed-Form Matrices for One-Dimensional 
Nodeless Variable Flux-Based Finite Element 


The closed-form expressions for the terms in the one-dimensional 
nodeless variable flux-based finite element matrices used in the thermal 
formulation, [M] and [D] are given by 

M(1 ,1) = M(2,2) = L/3 

M(1 ,2) = M(2,1) = L/6 

M(1 ,3) = M(2,3) = M (3,1) = M (3,2) = L/12 

M(3,3) = L/30 (A.1) 

D(1,1) = D(1,2) = -1/2 
D(2,1) = D(2,2) = 1/2 
D(3,1) = 1/6 

D(3,2) = -1/6 (A.2) 


For the structural formulation, the matrix [D 2 ] and the matrix [P] are given by 

D 2 (1.1) = D 2 (1.2) = -1/2 
D 2 (2,1) = D 2 (2,2) = 1/2 
D 2 (3,1) = 1/6 
D 2 (1,3) = -1/6 
D 2 (2,3) = 1/6 
D 2 (3,2) = -1/6 
D 2 (1,3) = -1/6 
D 2 (2,3) = 1/6 

D 2 (3,3) = 0 (A3) 
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P(1.1) = P(2,2) = E/L 
P(1,2) = P(2,1) = -E/L 

P(1,3) = P(2,3) = P(3,1) =P(3,2) = 0 
P(3,3) = E/3L 


(A.4) 
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Appendix B 

Closed-Form Matrices for Two-Dimensional 
Nodeless Variable Flux-Based Finite Element 

The determinant of the Jacobian matrix, |J|, is required for evaluation of 
the two-dimensional finite element matrices. It can be expressed in a simplified 

form as 


Ml- 


4 


iNiM Vj 
i=1 


(B.1) 


where Nj, i = 1 to 4 are the linear element interpolation functions. The terms v, 
are defined by 


V1 = ((x2-Xi)(y4-yi)-(x4-xi)(y2-yi))M 
V 2 = ((X2-Xi)(y3-y2)-(X3-X2)(y2-yi)V 4 
V3 = ((X3-X 2 )(y4-y3)-(X4-X3)(y3*y2))/4 

v 4 = ((X4-Xl)(y4-y3)-(X4-X3)(y4-yi))/ 4 ( B - 2 ) 

where xj and yj, i = 1 to 4 are the x- and y-coordinates of node i. The 
coefficients in the thermal finite element mass matrix are evaluated using the 
expression 


1 1 4 

M(i,j) = f fNjNj( £N k v k )dCdri (B.3) 

-1 -1 k=1 


where each coefficient is evaluated as 


M(1 ,1) = (9Vi + 3v 2 + v 3 + 3v 4 )/36 
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M(1 ,2) m (3vi + 3V2 + V3 + V4)/36 
M(1 ,3) = (vi + V2 + V3 + V4)/36 
M(1 ,4) = (3vi + V2 + V3 + 3v4)/36 
M(1 ,5) = (9v-| + 6v 2 + 2v 3 + 3v 4 )/180 
M(1 ,6) = (3vi + 3v 2 + 2v 3 + 2v4)/180 
M(1 ,7) = (3vi + 2v 2 + 2v 3 + 3v4)/1 80 
M(1 ,8) = (9v! + 3v 2 + 2v 3 + 6v 4 )/180 
M(2,1) = (3vi + 3 v 2 + V3 + V4)/36 
M(2,2) = (3vi + 9V2 + 3V3 + V4)/36 
M(2,3) = (vi + 3 v2 + 3v 3 + V4)/36 
M(2,4) = (vi + V2 + V3 + V4)/36 
M(2,5) = (6vi + 9v2 + 3v 3 + 2v4)/180 
M(2,6) = (3vi + 9v 2 + 6V3 + 2v4)/180 
M(2,7) = (2vi + 3v 2 + 3v 3 + 2v 4 )/180 
M(2,8) = (3vi + 3v2 + 2v 3 + 2v4)/180 
M ( 3 -1 ) = (vi + v 2 + v 3 + v 4 )/36 
M(3,2) = (vi + 3v2 + 3v 3 + V4)/36 
M(3,3) = (vi + 3 v 2 + 9V3 + 3v 4)/36 
M(3,4) = (vi + V2 + 3v 3 + 3v4)/36 
M(3,5) = (2 V1 + 3v 2 + 3v 3 + 2v 4 )/180 
M(3,6) = (2vi + 6V2 + 9V3 + 3v4)/180 
M(3,7) = (2v-| + 3v 2 + 9v 3 + 6v 4 )/180 
M(3,8) = (2vi + 2v 2 + 3v 3 + 3v 4 )/180 
M(4,1) = (3vi + v 2 + v 3 + 3v 4 )/36 
M(4,2) = (vi + V2 + V3 + V4)/36 
M(4,3) = (V! + v 2 + 3v 3 + 3v 4 )/36 
M(4,4) = (3vi + v 2 + 3v 3 + 9v 4 )/36 
M(4,5) = (3vi + 2v 2 + 2v 3 + 3v 4 )/180 
M(4,6) = (2v-| + 2v 2 + 3v 3 + 3v 4 )/180 
M(4,7) = (3v, + 2 v 2 + 6v 3 + 9v 4 )/180 
M(4,8) = (6vi + 2v 2 + 3V3 + 9v4)/180 
M(5,1) = (9vi + 6v 2 + 2V3 + 3v4)/180 
M(5,2) = (6vi + 9v 2 + 3v 3 + 2v 4 )/180 
M(5,3) = (2vi + 3v 2 + 3v 3 + 2v 4 )/180 
M(5,4) = (3vi + 2 v2 + 2V3 + 3v4)/180 
M(5,5) = (3vi + 3v 2 + V3 + V4)/180 
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M(5,6) = (6vi + 9v 2 + 6v 3 + 4v 4 )/900 
M(5,7) = (vi + v 2 + v 3 + v 4 )/1 80 
M(5,8) = (9vi + 6v 2 + 4v 3 + 6v 4 )/900 
M(6,1 ) = (3vi + 3v 2 + 2v 3 + 2v 4 )/180 
M(6,2) = (3vi + 9v 2 + 6v 3 + 2v 4 )/180 
M(6,3) = (2vi + 6v 2 + 9v 3 + 3v 4 )/180 
M(6,4) = (2vi + 2v 2 + 3v 3 + 3v 4 )/180 
M(6,5) = (6vi + 9v 2 + 6v 3 + 4v 4 )/900 
M(6,6) = (vi + 3v 2 + 3v 3 + v 4 )/180 
M(6,7) = (4vi + 6v 2 + 9v 3 + 6v 4 )/900 
M(6,8) = (vi + v 2 + v 3 + v 4 )/1 80 
M(7,1 ) = (3vi + 2v 2 + 2v 3 + 3v 4 )/180 
M(7,2) = (2vi + 3v 2 + 3v 3 + 2v 4 )/180 
M(7,3) = (2vi + 3v 2 + 9v 3 + 6v 4 )/180 
M(7,4) = (3vi + 2v 2 + 6v 3 + 9v 4 )/180 
M(7,5) = (vi + v 2 + v 3 + v 4 )/1 80 
M(7,6) = (4vi + 6v 2 + 9v 3 + 6v 4 )/900 
M(7,7) = (vi + v 2 + 3v 3 + 3v 4 )/1 80 
M(7,8) = (6vi + 4v 2 + 6v 3 + 9v 4 )/900 
M(8,1) = (9vi + 3v 2 + 2v 3 + 6v 4 )/180 
M(8,2) = (3vi + 3v 2 + 2v 3 + 2v 4 )/180 
M(8,3) = (2vi + 2v 2 + 3v 3 + 3v 4 )/180 
M(8,4) = (6vi + 2v 2 + 3v 3 + 9v 4 )/180 
M(8,5) = (9vi + 6v 2 + 4v 3 + 6v 4 )/900 
M(8,6) = (v 1 + v 2 + v 3 + v 4 )/1 80 
M(8,7) = (6vi + 4v 2 + 6v 3 + 9v 4 )/900 

M(8,8) = (3vi + v 2 + v 3 + 3v 4 )/180 (B.4) 


The coefficients in the thermal finite element matrices, [D x ] and [D y ], are 
evaluated using the expressions 

Dx(i.j) = [ -J, 2 ^)( N) T d^n (B.5) 
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1 1 


Dy(i ’ j) = .[ .[ (' J21 ^ + N (B.6) 


where each coefficient in D x (i,j) is evaluated as 


Dx(1.1) = (Y 2 - y4)/6 
Dx(1 ,2) = (2y2 - y3 - y4)/1 2 
Dx(1.3) = (y 2 - y 4 )/1 2 
D x(1 ,4) = (y 2 + y 3 - 2y 4 )/12 
D x(2,1) = (-2yi + y 3 + y 4 )/12 
D x (2,2) = (- yi + y 3 )/6 
d x( 2,3) = (- yi + 2y 3 - y 4 )/1 2 
D x (2,4) = (- yi +y3 )/l2 
Dx(3,1) = (-y 2 + y 4 )/12 
D x ( 3,2) = ( yi - 2y 2 + y4 )/1 2 
D x ( 3,3) = (- y2 + y 4 )/6 
D x (3,4) = (- yi - y2 + 2 y4 )/1 2 
D x (4.1) = (2 yi - y 2 - y 3 )/1 2 
D x (4,2) = ( yi - y3 )/1 2 
D x (4,3) = ( yi + y 2 - 2y 3 )/1 2 
D x (4,4) = ( yi - y3 )/6 
D x (5,1) = (-6 y i + 2y 2 + y 3 + 3y 4 )/72 
D x (5,2) = (-2 yi + 6y 2 - 3y 3 - y 4 )/72 
D x (5,3) = (- yi + 3y 2 - 2y 4 )[72 
D x (5,4) = (-3 yi + y2 + 2y 3 )/72 
D x (6,1) = (-3y 2 + y 3 + 2y 4 )/72 
D x (6,2) = (3 yi - 6y 2 + 2y 3 + y 4 )/72 
D x ( 6,3) = (- yi - 2y 2 + 6y 3 - 3y 4 )/72 
D x (6,4) = (-2 yi - y 2 + 3y 3 )/72 
D x (7,1) = (-2 y2 - y3 + 3 y4 )/72 
D x (7,2) = (2 yi - 3y 3 + y 4 )/72 
Dx(7,3) = ( yi + 3y 2 - 6y 3 + 2y 4 )/72 
D x (7,4) = (-3yi - y 2 - 2y 3 + 6y 4 )/72 
D x (8,1 ) = (6 yi - 3y 2 - y3 - 2y 4 )/72 
D x (8,2) = (3yi - 2y 3 - y 4 )/72 
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Dx(8,3) = (yi + 2y 2 - 3y 4 )/72 

D x (8,4) = (2yi + y 2 + 3y 3 - 6y 4 )/72 (B.7) 

and each coefficient in Dy(i,j) is evaluated as 

D y (1,1) = (*x 2 + x 4 )/6 
D y (1 ,2) = (-2x 2 + X 3 + x 4 )/1 2 
D y (1 ,3) = (-x 2 + x 4 )/12 
D y (1 ,4) = (*x 2 - x 3 + 2 x 4 )/12 
D y (2,1) = (2xi - x 3 - x 4 )/12 
D y (2,2) = (xi - x 3 )/6 
D y (2,3) = (xi -2x 3 + x 4 )/12 
D y (2,4) = (xi - x 3 )/1 2 
D y (3,1) = (x 2 - x 4 )/12 
D y (3,2) = (-xi + 2x 2 - x 4 )/1 2 
D y (3,3) = (x 2 - x 4 )/6 
D y (3,4) = (x-| + x 2 - 2x 4 )/1 2 
Dy(4,1 ) = (-2X1 + x 2 + x 3 )/1 2 
D y (4,2) = (-x-j + x 3 )/12 
D y (4,3) = (-x-i - x 2 + 2x 3 )/12 
D y (4,4) = (-xi + x 3 )/6 
D y (5,1 ) = (6xi - 2x 2 - x 3 - 3x 4 )/72 
D y (5,2) = (2xi - 6 x 2 + 3x 3 + x 4 )/72 
D y (5,3) = (xi - 3x 2 + 2x 4 )/72 
D y (5,4) = (3xi - x 2 - 2x 3 )/72 
D y (6, 1 ) = (3x 2 - x 3 - 2x 4 )/72 
D y (6,2) = (-3xi + 6 x 2 - 2x 3 - x 4 )/72 
D y (6,3) = (Xi + 2 x 2 - 6 x 3 + 3x 4 )/72 
D y (6,4) = (2xi + x 2 - 3x 3 )/72 
D y (7, 1 ) = (2x 2 + x 3 - 3x 4 )/72 
D y (7,2) = (-2xi + 3x 3 - x 4 )/72 
D y (7,3) = (-xi - 3x 2 + 6 x 3 - 2x 4 )/72 
D y (7,4) = (3xi + x 2 + 2x 3 - 6x 4 )/72 
D y (8, 1 ) = (-6X1 + 3x 2 + x 3 + 2x 4 )/72 
D y (8,2) = (-3xi + 2x 3 + x 4 )/72 
D y (8,3) = (-xi - 2x 2 + 3x 4 )/72 
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Dy(8,4) = (-2X1 - X 2 - 3x 3 + 6 x 4 )/72 


(B8) 


The boundary matrix, [B], for the thermal formulation is evaluated over the 
element surface (i.e., the element edge) using the one-dimensional 
interpolation functions. The coefficients in the thermal boundary matrix are 
given by 

B(1,1) = L/3 
B(1 ,2) = L/6 
B(2, 1 ) = L76 
B(2,2) = 173 
B(3,1 ) = 1712 

B(3,2) = 1712 (Bg) 

where L is the length of the element edge where the applied heating is defined. 

For the structural formulation, the coefficients in the finite element matrix, 
[Di] are given by 


Di(1,1) = (y 2 -y 4 )/6 
Di(1,2) = (2y 2 -y 3 -y 4 )/12 
Di (1 ,3) = (y 2 - y 4 )/1 2 
°i(1,4) = (y 2 + y 3 -2y 4 )/12 
Di(1,5) = 0 
Di( 1,6) = 0 
Di(1,7) = 0 
Di(1,8) = 0 
Di(1.9) = (-x 2 + x 4 )/6 
Di(1,10) = (-2x 2 + x 3 + x 4 )/12 
Di (1 ,1 1) = (-x 2 + x 4 )/12 
D-j (1 , 1 2) = (-x 2 - x 3 + 2x 4 )/12 
Di(2,1) = 0 

Di (2,2) = 0 
Di(2,3) = 0 
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Di (2,4) = 0 
Di(2,5) = (-x 2 + x 4 )/6 
Di( 2,6) = (-2x2 + x 3 + x 4 )/12 
Di (2,7) = (-x 2 + x 4 )/1 2 
Di (2,8) = (-x 2 - x 3 + 2 x 4 )/12 
Di( 2,9) = (y 2 * y 4 )/6 
Di(2,10) = (2y 2 -y 3 -y4)/12 
Di (2,1 1 ) « (y 2 - y 4 )/1 2 
Di (2,12) = (y 2 + y 3 - 2y 4 )/1 2 
D 1 (3,l) = (-2yi +y 3 + y 4 )/l2 
Di(3,2) = (-yi + y 3 )/6 
Di(3,3) = (-yi + 2y 3 - y 4 )/1 2 
Di(3,4) = (-yi + y 3 )/1 2 
Di (3,5) = 0 
Di (3,6) = 0 
Di (3,7) * 0 
Di (3,8) = 0 

Di (3.9) = (2xi - x 3 - x 4 )/12 
Di (3.1 0) = (xi - x 3 )/6 
D 1 (3,11) = (xi - 2x 3 + x 4 )/1 2 
Di(3,12) = (xi -x 3 )/12 
D-i (4,1 ) = 0 
Di(4,2) = 0 
Di(4,3) = 0 
Di(4,4) = 0 

D 1 (4,5) = (2xi -x 3 -x 4 )/12 
Di( 4.6) = (xi - x 3 )/6 
Di(4,7) = (xi-2x 3 + x 4 )/12 
D 1 (4.8) = (x 1 -x 3 )/12 
Di(4,9) = (-2yi +y 3 + y 4 )/12 
Di (4,1 0) = (-yi + y 3 )/6 
Di(4,1 1) = (-yi + 2y 3 - y 4 )/1 2 
Di(4,12) = (-yi +y 3 )/12 
Di(5,1) = (-y 2 + y 4 )/12 
Di(5,2) = (yi-2y 2 + y 4 )/12 
Di(5,3) = (*y 2 + y4)/6 
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Di (5,4) = (-yi - y 2 + 2y 4 )/12 

Di(5,5) = 0 

Di(5,6) = 0 

Di(5,7) = 0 

Di (5,8) = 0 

Di( 5,9) = (x 2 - x 4 )/12 

Di( 5,10) = (-x-j + 2x2 * x 4 )/12 

Di (5, 11) = (x 2 - x 4 )/6 

D-j (5,1 2) = (xi + x 2 - 2 x 4 )/1 2 

Di(6,1) = 0 

Di(6,2) = 0 

Di (6,3) = 0 

Di (6,4) = 0 

Di(6,5) = (x 2 - x 4 )/1 2 

Di(6,6)»(-xi +2 x 2 -x 4 )/12 

Di( 6,7) * (x 2 - x 4 )/6 

Di( 6,8) = (xi + x 2 - 2 x 4 )/1 2 

Di( 6,9) = (-y 2 + y 4 )/1 2 

Di(6,10) = (yi - 2y 2 + y 4 )/1 2 

Di(6,11) = (-y 2 + y 4 )/6 

Di(6,12) m (-yi - y 2 + 2y 4 )/12 

Di( 7.1) = (2yi - y 2 - y3)/1 2 

Di( 7,2) = (yi - y 3 )/1 2 

Di (7 ,3) = (yt + y 2 - 2y 3 )/1 2 

Di( 7,4) = (yi - y 3 )/6 

Di(7,5) = 0 

Di(7,6) = 0 

Di(7,7) = 0 

Di (7,8) = 0 

Di( 7,9) = (-2xi + x 2 + x 3 )/1 2 
Di(7,10) = (-xi + x 3 )/1 2 
Di(7,11) = (-xi -x 2 + 2x 3 )/12 
Di(7,12) = (-xi +x 3 )/6 
Di(8,1) = 0 

Di (8,2) = 0 

Di(8,3) = 0 
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Di (8,4) = 0 

Di(8,5) = (-2xi +x 2 + x 3 )/12 
Di(8,6) = (-xi +x 3 )/12 
Di(8,7) = (-xi -x 2 + 2x 3 )/12 
Di( 8,8) = (-xi + x 3 )/6 
Di(8,9) = (2yi-y 2 -y3)/12 
Di(8,l0) = (yi - y 3 )/l 2 
Di(8,11) = (yi + y 2 - 2y 3 )/1 2 
Di(8,12) = (yi-y 3 )/6 
Di (9.1 ) = (-6yi + 2y 2 + y 3 + 3y 4 )/72 
Di(9,2) = (-2yi + 6y 2 - 3y 3 - y 4 )/72 
Di(9,3) = (-yi + 3y 2 - 2y 4 )/72 
Di (9,4) = (-3yi + y 2 + 2y 3 )/72 
Di (9,5) = 0 
Dl(9,6) as 0 
Di(9,7) = 0 
Di(9,8) = 0 

Di(9,9) = (6xi - 2x 2 - x 3 - 3x 4 )/72 
Di (9,1 0) = (2xi - 6x 2 + 3x 3 + x 4 )/72 
Di(9,11) = (xi - 3x 2 + 2x 4 )/72 
Di (9,1 2) = (3xi - x 2 - 2x 3 )/72 
Di (10,1) = 0 

Di (10,2) = 0 
Di(10,3) = 0 
Di (10,4) = 0 

Di(10,5) = (6xi - 2x 2 - x 3 - 3x 4 )/72 
Di(10,6) = (2xi - 6x 2 + 3x 3 + x 4 )/72 
0^10,7) = (xi - 3x 2 + 2x 4 )/72 
Di(10,8) = (3xi - x 2 - 2x 3 )/72 
Di( 10,9) = (-6yi + 2y 2 + y 3 + 3y 4 )/72 
Di (1 0,1 0) = (-2yi + 6y 2 - 3y 3 - y 4 )/72 
Di (1 0, 1 1 ) = (-yi + 3y 2 - 2y 4 )/72 
Di(10,12) = (-3yi + y 2 + 2y 3 )/72 
Di (1 1,1) = (-3y 2 + y 3 + 2y 4 )/72 
Di (1 1 ,2) = (3yi - 6y 2 + 2y 3 + y 4 )/72 
Di (1 1 ,3) = (-yi - 2y 2 + 6y 3 - 3y 4 )/72 
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Di(11,4) = (-2yi -y 2 + 3y 3 )/72 

Di(11,5) = 0 

Di(11,6) = 0 

Di(11,7) = 0 

Di(11,8) = 0 

Di(1 1 ,9) = (3 x 2 - x 3 - 2 x 4 )/72 
Di( 1 1,10) = (-3xi + 6x 2 - 2x 3 - X4)/72 
Di(11,11) = (xi + 2x 2 - 6x 3 + 3x 4)/72 
D-|(1 1,12) = (2xi -»■ x 2 - 3 x 3 )/72 
Di (12,1 ) = 0 
Di (1 2,2) = 0 

Di(12,3) = 0 
Di (1 2,4) = 0 

Di (1 2,5) = (3X2 - x 3 - 2x 4 )/72 
Di (1 2,6) = (-3xi + 6x 2 - 2x 3 - x 4 )/72 
Di (1 2,7) = (xi + 2x 2 - 6x 3 + 3x 4 )/72 
Di(12,8) = (2xi + x 2 - 3x 3 )/72 
Di (1 2,9) = (-3y 2 + y 3 + 2y 4 )/72 
Di(12,10) = (3yi - 6y 2 + 2y 3 + y 4 )/72 
Di (1 2,1 1 ) = (-yi - 2y 2 + 6y 3 - 3y 4 )/72 
Di(12,12) = (-2yi -y 2 + 3y 3 )/72 
Di(13,1) = (-2y 2 - y 3 + 3y 4 )/72 
Di(13,2) = (2yi -3y 3 +y 4 )/72 
D-i(13,3) = (yi + 3y 2 - 6y 3 + 2y 4 )/72 
Di (1 3,4) = (-3yi - y 2 - 2y 3 + 6y 4 )/72 
Di(13,5) = 0 
Di (1 3,6) = 0 
Di(13,7) = 0 
Di(13,8) = 0 

Di(13,9) = (2x 2 + x 3 - 3x 4 )/72 
Di(13,10) = (-2xi + 3x 3 - x 4 )/72 
Di(13,1 1) = ( xi - 3x 2 + 6x 3 - 2x 4 )/72 
Di (1 3, 1 2) = (3xi + x 2 + 2x 3 - 6x 4 )/72 
Di(14,1) = 0 
Di (1 4,2) = 0 
Di (1 4,3) = 0 
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Di(14,4) = 0 

Di (1 4,5) = (2x 2 + x 3 - 3 x 4 )/72 
Di( 14,6) = (-2xi + 3x 3 - x 4 )/72 
Di( 14,7) = (-xi - 3 x 2 + 6x 3 - 2X4)172 
Di(14,8) = (3xi + x 2 + 2x 3 - 6 x 4)/72 
Di( 14,9) = (-2y 2 - Y3 + 3y 4 )/72 
Di (14,10) = (2yi -3y 3 + y 4 )/72 
Di (1 4,11) = (yi + 3y 2 - 6y 3 + 2y 4 )/72 
Di(14,12) = (-3yi - y 2 - 2y 3 + 6y 4 )/72 
Di(1 5,1 ) = (6yi - 3y 2 - y 3 - 2y 4 )/72 
Di(15,2) = (3yi -2y 3 -y 4 )/72 
Di (1 5,3) = (yi + 2y 2 - 3y 4 )/72 
Di(15,4) = (2yi + y 2 + 3y 3 - 6y 4 )/72 
Di (1 5,5) = 0 
Di(15,6) = 0 
Di (1 5,7) = 0 
Di(15,8) = 0 

Di(15,9) = (-6xi + 3x 2 + x 3 + 2x 4 )/72 
Di (1 5, 1 0) = (-3xi + 2x 3 + x 4 )/72 
Di( 15,1 1) = (-xi -2x 2 + 3x 4 )/72 
D! (15,1 2) = (-2X1 - x 2 - 3x 3 + 6x 4 )/72 
Di (16,1) = 0 

Di (16,2) = 0 
Di ( 1 6,3) * 0 
D t (16,4) = 0 

Di (1 6,5) = (-6xi + 3x 2 + x 3 + 2x 4 )/72 
Di(16,6) = (-3xi + 2x3 + x 4 )/72 
Di(16,7) = (-xi - 2x 2 + 3x 4 )/72 
Di(16,8) = (-2xi - x 2 - 3x 3 + 6x 4 )/72 
Di(16,9) = ( 6yi - 3y 2 - y3 - 2y 4 )/72 
Di(16,10) = (3yi -2y 3 -y 4 )/72 
Di (1 6,1 1 ) = (yi + 2y 2 - 3y 4 )/72 
Di (1 6,1 2) = (2yi + y 2 + 3y3 - 6y 4 )/72 


and the coefficients in the matrix [D2] are given by 


D 2 (1,1) = (y 2 - y4)/6 
□2(1,2) = (2y 2 - y 3 - y 4 )/1 2 
□ 2 (1.3) = (y 2 - y4)/12 
□ 2 (1 A) = (y 2 + y3 ■ 2y 4 )/l2 
□ 2 (1.5) = (4y 2 -y 3 -3y4)/72 
02(1,6) = (3y 2 - y 3 - 2y 4 )/72 
02(1,7) = (2y 2 + y 3 - 3y 4 )/72 
0 2 (1 ,8) = (3y 2 + y 3 - 4y 4 )/72 
0 2 (1 ,9) = 0 
D 2 (1,10) = 0 
O 2 (1,11) = 0 
0 2 (1 ,1 2) = 0 
0 2 (1 , 1 3) = 0 
D 2 (1,14) = 0 
0 2 (1 , 1 5) = 0 
0 2 (1 ,16) = 0 
D 2 (2,1) = 0 
D 2 (2,2) = 0 
D 2 (2,3) = 0 
02(2,4) = 0 
D 2 (2,5) = 0 
D 2 (2,6) = 0 
D 2 (2,7) = 0 
D 2 (2,8) = 0 
D 2 (2,9) = (-x 2 + x 4 )/6 
D 2 (2,10) = (-2 x 2 + x 3 + x 4 )/12 
D 2 (2, 1 1 ) = (-x 2 + x 4 )/1 2 
D 2 (2, 1 2) = (-x 2 - x 3 + 2 x 4 )/12 
D 2 (2,13) = (-4x 2 + x 3 + 3x 4 )/72 
D 2 (2,14) = (-3x 2 + x 3 + 2x 4 )/72 
02 ( 2 , 1 5) = (-2 x 2 - x 3 + 3x 4 )/72 
D 2 (2, 1 6) = (-3x 2 - x 3 + 4x 4 )/72 
D 2 (3,1) = (-2yi + y 3 + y 4 )/12 


f. 
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D 2 (3,2) = 

(-yi + 

ya)/6 


D 2 (3,3) = 

(*yi + 

2y 3 -y 4 )/i2 

D 2 (3,4) = 

(-yi + 

y3)/i 2 


D 2 (3,5) = 

(■4yi 

+ 3y 3 + 

y 4 )/72 

D 2 (3,6) = 

(-3yi 

+ 4y 3 - 

y 4 )/72 

D 2 (3,7) = 

(*2yi 

+ 3y 3 - 

y 4 )/72 

D 2 (3,8) = 

(*3yi 

+ 2y 3 + 

y 4 )/72 

D 2 (3,9)=( 

) 



D 2 (3,10) 

= 0 



D 2 (3,1 1 ) 

= 0 



D 2 (3,1 2) 

= 0 



D 2 (3,13) : 

= 0 



D 2 (3,14) 

= 0 



D 2 (3,15) 

= 0 



D 2 (3,16) 

= 0 



D 2 (4,1) = 

0 



D 2 (4,2) = 

0 



D 2 (4,3) = 

0 



D 2 (4,4) = 

0 



D 2 (4,5) = 

0 



D 2 (4,6) = 

0 



D 2 (4,7) = 

0 



D 2 (4,8) = 

0 



D 2 (4,9) = 

(2xi - 

x 3 - x 4 )/12 

D 2 (4,10) 

= (xi - 

x 3 )/6 


D 2 (4,1 1 ) 

= (xi - 

2x 3 + x 4 )/1 2 

D 2 (4,12) 

= (XI - 

x 3 )/12 


D 2 (4,13) 

= (4xi 

- 3x 3 - 

x 4 )/72 

D 2 (4,14) 

= (3xi 

- 4x 3 + 

x 4 )/72 

D 2 (4,15) 

= (2X! 

- 3x 3 + 

x 4 )/72 

D 2 (4,16) 

= (3x-| 

-2x 3 - 

x 4 )/72 

D 2 (5,1) = 

(-y 2 + 

y 4 )/12 


D 2 (5,2) = 

(yi - 2y 2 + y 4 )/12 

D 2 (5,3) = 

(-y 2 + 

y 4 )/6 


D 2 (5,4) = 

(-yi - 

y 2 + 2y 4 )/12 

D 2 (5,5) = 

(yi - 3y 2 + 2y 4 )/72 
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D 2 (5,6) = (yi - 4y 2 + 3y 4 )/72 

D 2 (5,7) = ( yi - 3y 2 + 4y 4 )/72 

D 2 (5,8) = (- yi - 2y 2 + 3y 4 )/72 

D 2 (5,9) = 0 

D 2 (5,10) = 0 

D 2 (5,1 1) = 0 

D 2 (5,12) =0 

D 2 (5,13) = 0 

D 2 (5,14) = 0 

D 2 (5,15) = 0 

D 2 (5,16) = 0 

D 2 (6,1) = 0 

D 2 (6,2) = 0 

D 2 (6,3) = 0 

D 2 (6,4) = 0 

D 2 (6,5) = 0 

D 2 (6,6) = 0 

D 2 (6,7) = 0 

D 2 (6,8) = 0 

D 2 (6,9) = (x 2 - x 4 )/1 2 

D 2 (6,10) = (-Xi + 2 x 2 - x 4 )/12 

D 2 (6, 11) = (x 2 - x 4 )/6 

D 2 (6, 1 2) = (xi + x 2 - 2 x 4 )/12 

D 2 (6,13) = (- Xl + 3 x 2 - 2x 4 )/72 

D 2 (6,14) = (-Xi + 4 x 2 - 3x 4 )/72 

D 2 (6,15) = (xi + 3 x 2 - 4x 4 )/72 

D 2 (6,16) = (xi + 2 x 2 - 3x 4 )/72 

D2(7,1) = (2yi - y 2 - y 3 )/1 2 

02(7,2) = (yi - y 3 )/1 2 

D 2 (7,3) = (yi + y 2 - 2y 3 )/12 

^2(7,4) = (yi - y 3 )/6 

D 2 (7,5) = (3yi - y 2 - 2y 3 )/72 

02(7,6) = (2yi + y 2 - 3y 3 )/72 

02(7,7) = (3yi + y 2 - 4y 3 )/72 

D 2 (7,8) = (4yi - y 2 - 3y 3 )/72 

D 2 (7,9) = 0 
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D 2 (7.10) 

= 0 


D 2 (7,11) 

= 0 


D 2 (7,12) 

= 0 


D 2 (7,13) 

= 0 


D 2 (7,14) 

= 0 


D 2 (7,15) 

= 0 


D 2 (7,16) 

= 0 


D 2 (8,1) = 

0 


D 2 (8,2) = 

0 


D 2 (8,3) = 

0 


D 2 (8.4) = 

0 


D 2 (8,5) = 

0 


D 2 (8,6) = 

0 


D 2 (8,7) = 

0 


D 2 (8,8) = 

0 


D 2 (8,9) = 

(-2xi 

+ X 2 + X3)/12 

D 2 (8,10) 

= (xi 

+ x 3 )/1 2 

D 2 (8,1 1) 

= (-xi 

- x 2 + 2 x 3 )/12 

D 2 (8,12) 

= (-Xi 

+ x 3 )/6 

D 2 (8,13) 

= (-3xi + x 2 + 2 x3)/72 

D 2 (8, 1 4) 

= (-2xi - x 2 + 3 x3)/72 

D 2 (8,15) 

= (-3xi - x 2 + 4 x3)/72 

D 2 (8, 1 6) 

= (-4xi + x 2 + 3 x3)/72 

D 2 (9,1) = 

(-6yi 

+ 2y 2 + y 3 + 3y4)/72 

D 2 (9,2) = 

(-2yi 

+ 6y 2 - 3y3 - y4)/72 

D 2 (9,3) = 

(-yi h 

- 3y2 - 2y4)/72 

D 2 (9,4) = 

(-3yi 

+ y 2 + 2y3)/72 

D 2 (9,5) = 

(-yi ^ 

- y 2 )/60 

D 2 (9,6) = 

(-yi ■» 

• 3y 2 - y 3 - y4)/1 44 

D 2 (9,7) = 

(-yi -» 

■ y 2 + y 3 - y 4 )/i 20 

D 2 (9,8) = 

(-3yi 

+ y 2 + y 3 + y4)/i44 

D 2 (9,9) = 

0 


D 2 (9,10) 

= 0 


D 2 (9,12) 

= 0 


D 2 (9,13) 

= 0 


D 2 (9,14) 

= 0 
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D 2 (9,15) = 0 
D 2 (9,16) = 0 
D 2 (10,1) = 0 
D 2 (10,2) = 0 
D 2 (10,3) = 0 
D 2 (10,4) = 0 
D 2 (10,5) = 0 
D 2 (10,6) = 0 
D 2 (1 0,7) = 0 

D 2 (10,9) = (6x-| - 2x 2 - x 3 - 3x 4 )/72 
D 2 (10,10) = (2xi - 6x 2 + 3x 3 + x 4 )/72 
D 2 (10,11) = (x 1 -3x 2 + 2x 4 )/72 
D 2 (10,12) = (3xi - x 2 - 2x 3 )/72 
D 2 (10,13) = (x 1 - x 2 )/60 
D 2 (10,14) = (xi - 3x 2 + x 3 + x 4 )/1 44 
D2(1 0, 1 5) = (xi - x 2 - x 3 + x 4 )/1 20 
D 2 (1 0,1 6) = (3xi - x 2 - x 3 - x 4 )/1 44 
D 2 (1 1 ,1 ) = (-3y 2 + y 3 + 2y 4 )/72 
D 2 (1 1 ,2) = (3yi - 6y 2 + 2y 3 + y 4 )/72 
D 2 (1 1 ,3) = (-yi - 2y 2 + 6y 3 - 3y 4 )/72 
D 2 (1 1 ,4) = (-2 yi - y 2 + 3y 3 )/72 
D 2 (1 1 ,5) = (yi - 3y 2 + y 3 + y 4 )/144 
D 2 (1 1,6) = (-y 2 + y 3 )/60 
D 2 (1 1 ,7) = (-yi - y 2 + 3y 3 - y 4 )/144 
D 2 (1 1,8) = (-yi - y 2 + y 3 + y 4 )/120 
D 2 (11,9) = 0 
D 2 (11,10) = 0 
D 2 (1 1,1 1) = o 
□ 2 (1 1,12) = 0 
D 2 (1 1,13) = 0 
D 2 (1 1,14) = 0 
D 2 (1 1,15) = 0 
D 2 (1 1 ,16) = 0 
D 2 (1 2, 1 ) = 0 
D 2 (12,2) = 0 
D 2 (12,3) = 0 
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D 2 (12,4) = 0 
D 2 (12,5) = 0 
D 2 (12,6) = 0 
D 2 (1 2,7) = 0 
D 2 (12,8) = 0 

D 2 (12,9) = (3x 2 - x 3 - 2x 4 )/72 

D 2 (1 2,10) = (-3X1 + 6 x 2 - 2x 3 - M)/72 

D 2 (1 2,11) = (xi + 2 x 2 - 6x3 + 3x 4 )/72 

D 2 (1 2,1 2) = (2xi + x 2 - 3x3)/72 

D 2 (12,13) = (-xi + 3 x 2 - x 3 - x 4 )/144 

D 2 (12,14) = (x 2 - x 3 )/60 

D 2 (12,15) = (xi + x 2 - 3x 3 + x 4 )/1 44 

D 2 (1 2,16) = (Xi + x 2 - x 3 - x 4 )/1 20 

D 2 (13,1) = (*2y 2 - y 3 + 3y 4 )/72 

D 2 (13,2) = (2yi -3y 3 +y 4 )/72 

D 2 (1 3,3) = (yi + 3y2 - 6y3 + 2y 4 )/72 

D 2 (1 3,4) = (*3yi - y 2 - 2y 3 + 6y 4 )/72 

d 2 (13,5) = (yi - y 2 - y3 + Y4)/120 

D 2 (13,6) = (yi + y2 * 3y 3 + y 4 )/144 

D 2 (13,7) = (*y3 + y 4 )/60 

D 2 (13,8) = (-yi - y2 - Y3 + 3y 4 )/1 44 

D 2 (1 3,9) = 0 

D 2 (13,10) = 0 

D 2 (13,1 1) = 0 

D 2 (13,12) = 0 

D 2 (13,13) = 0 

D 2 (1 3,1 4) = 0 

D 2 (1 3, 1 5) = 0 

D 2 (1 3,1 6) = 0 

D 2 (14,1) = 0 

D 2 (14,2) = 0 

D 2 (1 4,3) = 0 

D 2 (14,4) = 0 

D 2 (14,5) = 0 

D 2 (14,6) = 0 

D 2 (14,7) = 0 
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D 2 (14,8) = 0 

02(14,9) = (2x 2 + x 3 - 3x 4 )/72 
D 2 (14,10) = (-2xi + 3x 3 - x 4 )/72 
D 2 (1 4, 1 1 ) = (- Xl - 3x 2 + 6x 3 - 2x 4 )/72 
D 2 (14,12) = (3xi + x 2 + 2x 3 - 6x 4 )/72 
D 2 (14,13) = (-x-j + x 2 + x 3 - x 4 )/120 
D 2 (14,14) = (-xi - X2 + 3x3 - x 4 )/1 44 
D 2 (14,15) = (x 3 - x 4 )/60 
D 2 (14,16) = (xi + x 2 + x 3 - 3x 4 )/1 44 
0 2 (1 5,1 ) = (6yi - 3y 2 - y 3 - 2y 4 )/72 
D 2 (15,2) = (3yi - 2y 3 - y 4 )/72 
D 2 (15,3) = ^ + 2y 2 - 3y 4 )/72 
0 2 (1 5,4) = (2y t + y 2 + 3y 3 - 6y 4 )/72 
D 2 (15,5) = (3yi - y 2 - y 3 - y 4 )/144 
02(15,6) = (y-j + y 2 - y 3 - y 4 )/120 
0 2 (1 5,7) = ( yi + y 2 + y 3 - 3y 4 )/144 
0 2 (1 5,8) = (y! - y 4 )/60 
D 2 (15,9) * 0 
0 2 (1 5,1 0) = 0 
D 2 (15,1 1) = 0 
D 2 (1 5, 1 2) = 0 
0 2 (1 5, 1 3) = 0 
0 2 (1 5, 1 4) = 0 
D 2 ( 1 5, 1 5) = 0 
D 2 (15,16) = 0 
0 2 (1 6,1) = 0 
D 2 (16,2) = 0 
D 2 (16,3) = 0 
D 2 (16,4) = 0 
0 2 (1 6,5) = 0 
0 2 (1 6,6) = 0 
D 2 (16,7) = 0 
D 2 (16,8) = 0 

0 2 (1 6,9) = (-6X! + 3x 2 + X 3 + 2x 4 )/72 
02(16,10) = (-3X! + 2x 3 + x 4 )/72 * 

0 2 (1 6,11) = (- Xl - 2x 2 + 3 x 4 )/7 2 
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D 2 (16,12) = (-2X! - x 2 - 3x 3 + 6x4>/72 
D2(16,13) = (-3xi + x 2 + X3 + X4VI44 
D2(16,14) = (-xi - x 2 + X3 + X4)/120 
D2(16,15) = (-Xi - x 2 - X3 + 3x4)/ 144 

D2(16,16) = (-xi + X4V6O (B.11) 

As with the thermal formulation, the boundary matrix [B] used in the structural 
formulation is evaluated over the element surface using the one-dimensional 
element interpolation functions. The coefficients in the structural boundary 
matrix are given by 

B(1,1) = L/2 
B(1 ,2) = 0 
B(2,1) = 0 
B(2,2) = L/2 
B(3,1) = L/2 
B(3,2) = 0 
B(4,1) = 0 
B(4,2) = U2 
B(5,1) = L/6 
B(5,2) = 0 
B(6,1) as 0 

B(6,2) st l_/6 (B.12) 

where L is the length of the element edge where the applied pressure is 
defined. The mechanical stress component vector { 01 } is defined in equation 

(4.50) as { ci } = [P]{5>, where the coefficients in the matrix [P] are given by 

P( 1, 1) = cn(y2-y4)/down1 
P( 1 , 2) = Ci2(-x 2 + X4)/down1 
P( 1 , 3) = cn(-yi + y 4 )/down1 
P( 1, 4) = -Ci2(-xi + X4)/down1 
P( 1. 5) = 0 
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P( 1, 6) = 0 

P( 1, 7) = *cn(-yi + y2)/down1 
P( 1 , 8) = Ci2(-xi + X2)/down1 
P( 1, 9) = Cn(-yi + y4)/down1 
P( 1,10) = -Ci2(-xi + X4)/down1 
P(1,11) = 0 
P( 1,12) = 0 
P( 1,13) = 0 
P( 1 ,14) = 0 

P( 1,15)= -Cn(-yi + y2)/down1 
P( 1,16) = Ci 2 (-xi + X2)/down1 
P( 2, 1) = -cn(-y 2 + y3)/down2 
P( 2, 2) = Ci 2 (-X 2 + X3)/down2 
P( 2, 3) = cii(-yi + y3)/down2 
P( 2, 4) = Ci 2 (xi - X3)/down2 
P( 2, 5) = -ci i(-yi + y2)/down2 
P( 2, 6) = Ci 2 (-xi + X2)/down2 
P(2, 7) = 0 
P( 2, 8) = 0 

P( 2, 9) = -ci i(-y 2 + y3)/down2 
P( 2,10) = Ci2(-X2 + X3)/down2 
P( 2,1 1) = -cii(-yi + y2)/down2 
P( 2,12) = Ci2(-xi + X2)/down2 
P( 2,13) = 0 
P( 2,14) = 0 
P( 2,15) = 0 
P( 2,16) = 0 
P( 3, 1 ) = 0 
P( 3, 2) = 0 

P( 3, 3) = cn(y 3 - y4)/down3 
P( 3, 4) = -Ci 2 (x 3 - X4)/down3 
P( 3, 5) = cn(-y 2 + y4)/down3 
P( 3, 6) = Ci 2 (x 2 - X4)/down3 
P( 3, 7) = -cn(-y 2 + y3)/down3 
P( 3, 8) = Ci2(-X2 + X3)/down3 
P( 3, 9) = 0 
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P( 3,10) = 0 

P( 3,1 1) = cn(y 3 - y 4 )/down3 
P( 3,12) = -ci2(x 3 - x 4 )/down3 
P( 3,13) = -cn(-y2 + y 3 )/down3 
P( 3,14) = ci2(-X2 + x 3 )/down3 
P( 3,15) = 0 
P( 3,16) = 0 

P( 4, 1) = cn(y 3 - y 4 )/down4 
P( 4, 2)=-ci 2 (x 3 - x 4 )/down4 
P( 4. 3) = 0 
P( 4, 4) = 0 

P( 4, 5) = cn(-yi + y 4 )/down4 
P( 4 , 6) = -ci2(*xi + x 4 )/down4 
P( 4, 7) = cn(yi -y 3 )/down4 
P( 4, 8)= ci 2 (-xi + x 3 )/down4 
P( 4, 9) = 0 
P( 4,10) = 0 
P( 4,1 1) = 0 
P( 4,12) = 0 

P( 4,13) = cn(-yi + y 4 )/down4 
P( 4,14) = -ci 2 (-xi + x 4 )/down4 
P( 4,15) = cn(y 3 - y 4 )/down4 
P( 4,16) = -ci 2 (x 3 - x 4 )/down4 
P( 5, 1) = C2i(y2 - y 4 )/down1 
P( 5, 2) = C2i(-yi + y 4 )/down1 
P( 5, 4) = -c 2 2 (-xi + x 4 )/down1 
P( 5. 5) = 0 
P( 5, 6) = 0 

P( 5, 7) = -c 2 i(-yi + y2)/downl 
P( 5, 8) = c 2 2(-xi + x 2 )/down1 
P( 5, 9) = c 2 i(-yi + y 4 )/downl 
P( 5,10) = -c 2 2(-xi + x 4 )/down1 
P( 5,11) = 0 
P( 5,12) = 0 
P( 5,13) = 0 
P( 5,14) = 0 
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P( 5,15) = -C 2 i(-yi + y2)/down1 
P( 5,16) = C 22 (-xi + X2)/down1 
P( 6, 1) = -C 2 i(-y 2 + y3)/down2 
P( 6, 2) = C 22 (-X 2 + X3)/down2 
P( 6, 3) = C 2 i(-yi + y3)/down2 
P( 6, 4) = C 22 (xi - x 3 )/down2 
P( 6, 5 ) = -C 2 i(-yi + y 2 )/down 2 
P( 6, 6) = C 22 (-xi + X2)/down2 
P( 6, 7) = 0 
P( 6, 8) = 0 

p ( 6, 9) = -C 2 i(*y 2 + y3)/down2 
P( 6,10) = C 22 (-X 2 + X3)/down2 
P( 6,1 1) = -C 2 i(-yi + y2)/down2 
P( 6,12) = C 22 (-xi + X2)/down2 
P( 6,13) = 0 
P( 6,14) = 0 
P( 6,15) = 0 
P( 6,16) = 0 
P( 7, 1 ) = 0 
P( 7, 2) = 0 

P( 7, 4) = -C 22 (X 3 - X4)/down3 
P( 7, 5) = C 2 i(-y 2 + y4)/down3 
P( 7, 6) = C 22 (X 2 - x 4 )/down3 
P( 7, 7) = -C 2 i(-y 2 + y3)/down3 
P( 7, 8) = C 22 (-X 2 + X3)/down3 
P( 7, 9) = 0 
P( 7,10) = 0 

P( 7, 1 1 ) = C 21 (y 3 - y4)/down3 
P( 7,12) = -C 22 (X 3 - X4)/down3 
P( 7,13) = -C 2 i(-y 2 + y3)/down3 
P( 7,14) s C 22 (-X 2 + X3)/down3 
P( 7,15) = 0 
P( 7,16) = 0 

P( 8, 1) = C 2 i(y 3 - y4)/down4 
P( 8, 2) = -C 22 (x 3 - X4)/down4 
P( 8, 3) = 0 
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P( 8, 4) = 0 

P( 8, 5) = C 2 i(*yi + y4)/down4 
P( 8, 6) = -C 22 (-xi + X4)/down4 
P( 8, 7) = C 2 i(yi - y3)/down4 
P( 8, 8) = C 22 (*xi + X3)/down4 
P( 8, 9) = 0 
P( 8,10) = 0 
P( 8,1 1) =0 
P( 8,12) = 0 

P( 8,13) = C 2 i(-yi + y4)/down4 
P( 8,14) = -C 22 (-xi + X4)/down4 
P( 8,1 5) = C 21 (y 3 - y4)/down4 
P( 8,16) = -C 22 (X 3 - X4)/down4 
P( 9, 1) = C 33 (-X 2 + X4)/down1 
P( 9, 2) = C 33 (y 2 - y4)/down1 
P( 9, 3) = -C 33 (*xi + X4)/down1 
P( 9, 4) = C 33 (-yi + y4)/down1 
P( 9, 5) = 0 
P( 9, 6) = 0 

P( 9, 7) = C 33 (-x-| + X2)/down1 
P( 9, 8) = -C 33 (-yi + y2)/down1 
P( 9, 9) = -C 33 (-xi + x4)/down1 
P( 9,10) = C 33 (-yi + y4)/down1 
P( 9,1 1) = 0 
P( 9,12) = 0 
P( 9,13) = 0 
P( 9,14) = 0 

P( 9,15) = C 33 (-xi + X2)/down1 
P( 9,16) = -C33(-yi + y 2 )/down 1 
P(10, 1) = C 33 (-x 2 + X3)/down2 
P(10, 2) = -C 33 (-y 2 + y3)/down2 
P(10, 3) = C 33 (xi - x 3 )/down2 
P(10, 4) = C 33 (-yi + y3)/down2 
P(10, 5) = C 33 (-xi + X2)/down2 
P(10, 6) = -C 33 (-yi + y2)/down2 
P(10, 7) = 0 
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P(10, 8) = 0 

P(10, 9) = C 33 (-X 2 + x 3 )/down2 

P(1 0, 1 0) = -C 33 (-y 2 + y3)/down2 

P(1 0,1 1 ) = C 33 (-xi + X2)/down2 

P(10,12) = -c 33 (-yi + y2)/down2 

P(10,13) = 0 

P(1 0, 1 4) = 0 

P(1 0,1 5) = 0 

P(1 0, 1 6) = 0 

P(11.1) = 0 

P(1 1, 2) = 0 

P(1 1, 3) = -C 33 (X 3 - x 4 )/down3 
P(1 1 , 4) = C 33 (y 3 - y 4 )/down3 
P(1 1 , 5) = C 33 (X 2 - x 4 )/down3 
P(1 1 , 6) = C 33 (-y 2 + y 4 )/down3 
P(1 1, 7) = C33(-X2 + X3)/dOWn3 
P(1 1.8) = -C 33 (*y 2 + y3)/down3 
P(1 1 , 9) = 0 
P(1 1,10) = 0 

P(1 1,11) = -C 33 (x 3 - x 4 )/down3 
P(1 1,12) = C 33 (y 3 - y 4 )/down3 
P(1 1,13) = C 33 (-X 2 + x 3 )/down3 
P(1 1,14) = -C33(-V2 + y3)/down3 
P(1 1,15) = 0 
P(1 1,16) = 0 

P(12, 1) = -c 33 (x 3 - x 4 )/down4 
P(12, 2) = c 33 (y 3 - y 4 )/down4 
P(12, 3) = 0 
P(12, 4) = 0 

P(12, 5) = -c 33 (-xi + x 4 )/down4 
P(12, 6) = c 33 (-yi + y 4 )/down4 
P(12, 7) = C 33 (-x^ + x 3 )/down4 
P(12, 8) = c 33 (yi - y 3 )/down4 
P(12, 9) = 0 
P(1 2, 1 0) = 0 
P( 1 2, 1 1 ) = 0 
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P(1 2,1 2) = 0 

P(12,13) = -C33(-Xi + x 4 )/down4 
P(1 2,1 4) = C33(-yi + y 4 )/down4 
P(1 2,1 5) = -C33(X3 * x 4 )/down4 

P(1 2,1 6) = C33(Y3 - y 4 )/down4 (B.13) 

where downl, down2, down3, and down4 are defined by 

downl = (x 2 *xi)y 4 +(xi-x 4 )y 2 +(x 4 -x 2 )yi 
down2 = (x2-xi)y3+(xi-X3)y2+(X3-X2)yi 
down3 = (X3-X2)y 4 +(X2-X 4 )y3+( X 4 -X3)y2 

down4 = (x 3 -xi)y 4 +(xi-x 4 )y 3 +(x 4 -x 3 )yi (B.14) 

and cjj, where i = 1 to 3 and j = 1 to 3, are the coefficients in the material elastic 
constant matrix defined in equation (4.40) for plane stress problems and 
equation (4.41) for plane strain problems. The coefficients in the thermal stress 
component vector { 02 } are given by 


o 2 (1)= a(cn +C12KT1 -T 0 ) 
02(2)= a (cii + C12MT2 -T 0 ) 
<*2 (3) = a (cn + C12XT3 - T 0 ) 
02 (4) = a (cn + ci2)(T 4 - T 0 ) 
02 (5) = a (cn + C12HT5) 

02 (6) = a (cn + Cl2)(T6) 

02 ( 7 ) = a (C11 +ci 2 )(T 7 ) 

02 (8) = a (C11 + C12XT8 ) 

02 ( 9 ) = a (C21 + C22HT1 - T 0 ) 
02 ( 10 ) = a (C21 + C22HT2 - To) 
02 ( 11 ) = a (C21 + C22HT3 - To) 
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02 (12) = a (C21 + C22KT4 - To) 

02 (13) = a (c 2 i + c 22 )(T 5 ) 

°2 (14) = a (C21 + C22XT6) 

02 (15)= a(c 2 i +C22KT7) 

o 2 (16)= a (c 2 i + c 22 )(T 8 ) (B15) 

where cij, 1 = 1 to 3 and j = 1 to 3, are the coefficients in the material elastic 

constant matrix defined in equation (4.40) for plane stress problems and 
equation (4.41) for plane strainjjroblems. For plane stress problems, « 
and for plane strain problems a = a(1+v). 
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